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ABSTRACT 


In this thesis the performance of a step quantum well infrared photodetector, designed 
by Kevin Lantz (June 2002) and experimentally studied by Michael Touse (September 2003) 
and Yeo Hwee Tiong (December 2004), was simulated in Matlab using the transfer matrix 
method. The results, obtained by the Matlab program, are compared with the experimental 
results, in an attempt to make inferences about the optimum way of designing QWIP detec- 


tors. 


Simulation of the above implies numerical solution of the Schrédinger equation, us- 
ing algorithms and methods, which give accurate results. In our approach, the transfer matrix 
method (TMM) was used with exponentials and Airy functions to represent the solutions to 
Schrédinger equation under zero and non-zero bias, respectively. The calculated results were 
compared with the experimental data and found to provide a good agreement which validated 


the accuracy of the model employed. 


In the final section of the thesis we examine and simulate in Matlab the application of 
the extended Kalman filtering (EKF) to an infrared photodetector as a target tracking mecha- 
nism to both maneuvering and non-maneuvering targets. When we used one sensor for track- 
ing, the results were reliable provided that the target did not maneuver. In the case of a ma- 
neuvering target the results were significantly improved when we used both sensors for 


tracking. 
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I. INTRODUCTION 


A. QUANTUM WELL INFRARED PHOTODETECTORS (QWIPS)-BRIEF 
REVIEW 


In the recent years there has been a considerable increase in research activities 
towards the development of quantum well based infrared sensors [1-4]. The use of quan- 
tum wells made of larger bandgap materials, can overcome the difficulties associated 
with smaller bandgap materials, typically used in conventional infrared detectors [1]. 
The smaller bandgap materials are relatively unstable, and difficult to process, resulting 
in non-uniformities when focal plane arrays are made [1]. The infrared detectors have 


many commercial and military applications. 


The operation of quantum well infrared photodetectors (QWIPs) is based on the 
excitation of bound electrons in a quantum well by infrared, as illustrated in Figure 1.1a. 
The quantum wells are formed by sandwiching a small bandgap material between large 
bandgap materials, as shown in Figure 1.1b. It is possible to use the quantum well 
formed either in the conduction band or valence band to fabricate QWIPs. In the case of 
conventional infrared detectors, the electrons in the valence band are excited across the 
bandgap, while in QWIP detectors it is necessary to dope the quantum well to populate 
the ground state using either n-type (for the conduction band well), or p-type (for the va- 


lence band well). 





(b) 
Ey 
Figure 1.1 (a) excitation of bound electrons due to incident photon energy ia . (b) 
Quantum well formation by sandwiching the appropriate bandgap material (From Ref. 


[5].) 


In the design of quantum wells, the parameters, such as its dimension, and the 
composition of each structure material used, allow us to manipulate the behavior of the 
photoexcited carriers, so that they escape from the potential wells (quantum leak) and are 


collected as current by the application of an external bias, as shown in Figure 1.2. 





photocurrent 





“dark current” 
mechanisms 





position 


Figure 1.2 Conduction band diagram of a QIWP in an external electric field. Excited 
electrons escape from the ground state creating a photocurrent (From Ref. [4].) 


The layers of QWIP detectors are primarily grown by molecular beam epitaxy 
(MBE), which gives thickness control in atomic scale [2, 4]. For achieving good material 
quality, the lattice constant of different materials used in QWIPs should be nearly the 


same. This avoids the dislocations generated by the lattice mismatch. 


The usual wavelength coverage of QWIPs is between 3 and 20 microns [1, 8]. By 
controlling the barrier height and well width, it is possible to adjust the energy level sepa- 
ration and, hence, the wavelength dependence of the detector response. For example, the 


longer wavelength response can be achieved using shallow quantum wells. 


QWIPs can be categorized according to the electron resulting state in four types: 
bound to bound (B-B), bound to quasibound (B-QB), bound to continuum (B-C), and 


bound to miniband (B-M), as schematically shown in Figure 1.3. 
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Figure 1.3 Energy band diagram showing B-B, B-QB, B-C, and B-M QWIP struc- 
tures (From Ref. [4].) 


The main factor which limits the QWIP performance is the dark current (the cur- 
rent that flows through a biased detector in the dark with no photons impinging on it) , 
which plagues all infrared detectors. At temperatures above 45 K, the dark current of the 
QWIP is entirely dominated by classic thermionic emission of ground state electrons di- 
rectly out of the well into the energy continuum. There are several theoretical techniques 
of minimizing it [3], most importantly the positioning of the first excited state to align 


with the barrier. 


Recent studies indicate that multicolor operation of QWIPs is possible by simply 
stacking different types of quantum wells. A four-color QWIP was reported [5], based on 
stacked InGaAs / AlGaAs andGaAs/ AlGaAs multi-quantum wells. In addition, a great 
deal of QWIP focal plane arrays have been analyzed and demonstrated by many groups 
[6-8], showing excellent performance in LWIR atmospheric window. 

B. PURPOSE OF THIS THESIS 

The purpose of this thesis was to simulate in Matlab the performance of symmet- 
ric and asymmetric quantum well structures under both zero bias and non-zero bias. It 
uses a design made in earlier thesis by Kevin R. Lantz [25] and experimental measure- 


ments made by Michael P. Touse [26] and Yeo Hwee Tiong [27]. The results obtained by 
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the MATLAB program are compared with the experimental results, in an attempt to make 
inferences about the optimum way of designing QWIP detectors. Simulation of the above 
implies numerical solution of the Schrédinger equation, using algorithms and methods 
which give accurate results. Various methods have been presented to solve the 
Schrédinger equation numerically including the WKB approximation, the Monte Carlo 
method and the finite element method, which all basically deal with solving second order 
differential equations. In our approach, the transfer matrix method will be used [12, 16] 
with exponential and Airy functions to represent the solutions to Schrédinger equation 


under zero and non-zero bias, respectively. 


In the final section of the thesis, we will examine and simulate in Matlab the ap- 
plication of the extended Kalman filtering (EKF) to an infrared photodetector as a target 
tracking mechanism, to both maneuvering and non-maneuvering targets. Using a recur- 
sive type algorithm such as EKF, we will examine the bearing only tracking (BOT) prob- 
lem classically. The main drawback of BOT is that the research is closed, meaning that 
the publication available is limited. 

C. MILITARY RELEVANCE 

QWIPs have a great deal of military and security applications. The most common 
are night vision, remote sensing and satellite IR imaging, FLIR, surveillance, targeting in 
a variety of terrains (air sea ground) and weapon delivery [4,7]. They are also widely 
used in search and rescue situations, mine detection, missile seeker formulation, smart 
munitions, weapon sights, preventive maintenance, non-destructive testing, medical im- 
aging and surveillance [4]. Finally, the bearing-only-tracking application examined in 
the last chapter of the thesis, is considered a modern real-world problem featuring the ad- 


vantage of passive target acquisition and surveillance. Since modern military needs are 


oriented in systems which can track targets without being noticed, the use of a modern 
QWIP sensor combined with a classic tracking algorithm, such as EKF, could result to a 
reliable surveillance system. 
D. THESIS OUTLINE 

The present thesis starts with a general discussion of the necessity of QWIPs, fol- 
lowed by a brief review of how these devices are formulated and where they find applica- 


tions of military interest. Next, the reader is led, through the mathematical background, to 
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the simulation program, and the comparison to the laboratory case. Finally a “real world” 


application is investigated through BOT. 
In particular, the outline of the present thesis is as follows: 


Chapter II discusses the Schrédinger equation and the formulation of the transfer 
matrix method. The Schrédinger equation is applied in solving the unbiased and biased 
quantum well, and the transfer matrix method is applied in order to expand the solution in 


more than one region of the structure. 


Chapter II provides information about the constructed Matlab code, and the labo- 
ratory semiconductor, while simulations are presented to show the QWIP behavior to 
several conditions of bias. Ideas such as probability current and bound-to-continuum tran- 


sitions are presented, and finally the comparison to the laboratory case takes place. 


Chapter IV is focused into the BOT problem. The used initialization of the state 
vector is presented, as long as the EKF algorithm, with emphasis into bearing tracking. 


Finally, we present and discuss several cases of sensor-target geometry. 


Chapter V is a summary of the results, followed by suggestions for probable fu- 


ture research. 
Appendix A includes the Matlab coding used to model the QWIP. 
Appendix B includes the Matlab coding used to model the BOT problem. 


Appendix C is an overview of Airy functions and their asymptotic expressions 


used to model the QWIP in the low bias cases. 


Finally, Appendix D includes the complete set of resulting graphs of the BOT 


problem. 
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I. MATHEMATICAL AND PHYSICAL BACKGROUND 


A. INTRODUCTION 

The infrared photodetectors examined in this thesis are quantum mechanical de- 
vices. Modeling this type of device requires solution of the key quantum mechanical 
equation, Schrédinger’s equation. Given the potential and boundary conditions, we de- 
termine the shape and the temporal evolvement of the wave function by solving this sec- 
ond-order, partial differential equation. Due to the one dimensional nature of quantum 
well structures used in QWIPs, it is sufficient to solve the Schrédinger equation along the 


growth direction of the layers. 


Over the years, various methods have been used for numerical solution of the 
Schrédinger equation, namely, the WKB approximation [9], the variational calculation 
method [10, 11], the Monte Carlo method [12], the finite element method (FEM) [13], 
and the transfer matrix method (TMM) [14]. In the present work, we employed the trans- 
fer matrix method (TMM) due to its flexibility in analyzing complex quantum well struc- 
tures under external bias. 

B. SCHRODINGER EQUATION 
The one dimensional time-dependent Schrédinger equation for an electron can be 


written as: 


2 
wan suy Yat) = ih W(,1), (2.1) 
2m, Ox ot 


where fi is the reduced Plank’s constant, m, is the electron mass, x is the coordinate 


(usually taken along the direction of the growth), ‘¥ is the wavefunction and U is the po- 


tential energy. 


Since the quantum well potential is independent of time, it is possible to separate 
the spatial and temporal dependencies with the spatial dependence satisfying the time- 
independent Schrédinger equation: 

ho 


~— = P(x) +U (x) P(x) = EV(4). (2.2) 
2m Ox 


C. EFFECTIVE MASS MODEL 
The periodic potential of the atoms in semiconductor materials influences the pro- 
pagation of an electron inside. In order to include the effect of the periodic potential, the 
electron mass m is usually replaced by an effective mass m’ , and treats the electron as a 
free particle [1]. The effective mass m is related to the band structure by: 
* hr 
OPE dR 


where fiis the Plank’s constant, FE is the conduction band energy and k the wavevector. 


(2.3) 


The QWIPs are usually made using potential wells in the conduction band well 
due to smaller effective mass which gives higher absorption [1]. For the materials used 
in this thesis (InGaAs , AlGaAs ) the effective mass of electrons for various composi- 


tions is shown in Figure 2.1 








Effective mass in Kgr 














0.5 
Doping percentage x/100 


Figure 2.1: Effective mass as a function ofx for Al.Ga,_.As and In,Ga,_,As. 


D. TRANSFER MATRIX METHOD (TMM) 

1. General 

The transfer matrix method is a relatively simple approach for solving Schrédin- 
ger equation for quantum well structures. It is generally used for a number of problems 
mostly involving second-order, partial differential equations (e.g., propagation character- 


istics in planar waveguides, spontaneous emission in DFB lasers) [12-14, 16]. 


Comparing this method with the other methods usually used for solving the same 
problems, TMM has the advantage both of simplicity and usage in low processing ability 
computers. Its simplicity mainly involves multiplication of 2x2 matrices which are 
straightforward to implement on a computer. We will examine two cases using the TMM, 
quantum well structures with and without the application of an external electric field. 

2: Quantum Well Structure without Applied Electric Field 

We assume an arbitrary multilayered quantum well structure as shown in Figure 


Did 





Figure 2.2: | Multilayered quantum well structure with zero applied bias. 


The time independent Schrédinger equation for each square quantum well region 
of the multilayered quantum well structure reads: 
he OW (x) 


~=——— + VV," (a)¥ (x) = EY (x). (2.4) 
2m, Ox, . ‘ 





where /7is the reduced Plank’s constant, m, is the electron effective mass in the region j, 


xis the coordinate (usually taken along the direction of the growth), ‘V is the wave- func- 


tion, V;is the potential energy in region j, and E stands for the associated electron en- 


ergy. The general solution of Equation (2.17) in each region is a superposition of the left 


and the right traveling waves and is given by: 

Gj Ce" iDe (2.5) 
where k, is implemented from the separation of variables and represents the wavenumber 
and C,,D,are arbitrary constants. The value of k, depends on the region effective mass 


m, and the potential height V, in the region j as: 


2m;(E-V;) 


k.= 


: : (2.6) 


At each boundary the probability current should be continuous which implies: 





¥ (x)= Fal), (2.7) 
ld l d 
m, dx pee (x; oe - —— | ¥ a), (2.8) 


where xis the coordinate at the boundary between j and j +1 layers. 


Applying the solution of Equation (2.4), given by Equation (2.5), to Equations 
(2.7) and (2.8), we can find after some algebra: 








_ eat + esa enitie Ke mC, + Kary etn Ke a nD (2.9) 
De 2k 
jl; Lam. 
k..m.—k.m. k. m, “+km. 
1 1 ik j ik 1 j+1 —ik jy +ik 4x 
i epee See a 7 fe ad Ge i a 7 cee ab atc De (2.10) 
m, 


jt jt 


The results in Equations (2.9) and (2.10) can be expressed in a matrix form as 


Cia C; 
=M | ‘|, (2.11) 
D.. a. 





where M is given by as 


1+ i Ps Ki) X; t= i pays 
1 jl; ll; 


M = ane (2.12) 
j 2 P * : : 
a lia oiitkwdsy | 4 Pea diet oka 
kK jm; K ym, 


Repeatedly applying Equation (2.11), it is possible to relate the coefficients if the 


outermost layers as: 


Cy C, m, my, \\C, 
=M,M,, one M,M, = . (2.13) 
Dy D, My, My }| D, 


In the case of bound states, the wavefunction must decrease exponentially to zero 
at the outermost boundaries. This implies that exponential growing coefficients in Equa- 


tion (2.5), namely (D,,,C, ), should equate to zero. Under this condition Equation (2.13) 


eee 
= . (2.14) 
0 My, My, )| D, 


It can be easily seen from Equation (2.14) that, to satisfy the boundary conditions, 


reduces to: 


we demand the matrix element m,, = 0. Since m,, depends only on energy FE, the bound 
states can be found by plotting m,, (E ) and identifying the zero crossings. 
3: Quantum Well Structure under Applied Electric Field 


a. General 

During the QWIP operation, it is necessary to apply external bias to extract 
photoexcited electrons. Under an external electric field, the quantum well potential tilts, 
as illustrated in Figure 2.3. The amount of the tilt depends on the strength and the direc- 
tion of the applied field. The bias alters the energy levels in the quantum well structure, 
and hence, the photoresponse. In the section following, the effect of an external electric 
field on QWIP operation is analyzed by solving Schrédinger equation with a linear poten- 


tial in addition to the quantum well potential. 
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Figure 2.3 Multilayered quantum well structure with forward applied bias (left), and 
reverse applied bias (right). 


b. Solution of the Schrodinger Equation in a Linear Potential 


We first consider the solution of Schrédinger equation in a linear potential 


as shown in Figure 2.4. 





Figure 2.4: _ Linear potential. 


Using Schrédinger’s Equation (2.2) and substituting the value of 
U(x) =—ax we can get to the following form: 
i? Oo 


—~——'P(x)-ax¥(x) = EV (x). (2.15) 
2m Ox 


It is possible to convert Equation (2.15) to a standard differential equation 


by using the following change of variables. 


ae (2.16) 


x= Byty => y=— 
B 
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where @ and y have units of length, while y is considered dimensionless. 


Differentiating (2.16) with respect to x and substituting back into (2.15) 














yields: 
@ | 4-7 |) dyad) ey) 1d (2.17) 
dx\ B dx dy\ B B dy” 
1 d’ 2ma 2m 
2, 2 3 2 
. w+ ue yer IE (ay +E)¥ =0. (2.19) 


The variables @ andy can be chosen as given in Equation (2.16) to further 


simplify Equation (2.19), so: 








ie Sy, 
a 
2map° he ie ee) 
=-l>f=- : 

he e 2ma 
Using the parameters as defined in Equation (2.20), the Equation (2.19) 

can be transformed into the Airy’s differential equation: 

d’ 
—~¥-y¥=0, (2.21) 


which has two linearly independent solutions and are given in terms of Airy functions of 
the first and the second kind ( Ai and Bi). The general solution can be written as a linear 
combination of the two Airy function multiplied by two arbitrary constants C and D [16]: 
Y(y) = CAi(y)+DBi(y). (2.22) 
Figures (2.5) and (2.6) show the Airy functions of the first kind, Ai, and 


the second kind, Bi, as a function of y for both positive and negative directions. 
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Airy functions of the second kind. 
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Figure 2.6: 


In Equation (2.16) by replacing the values of yand / we get: 


E 
x+— 
a 


i 1/3 
2ma 
where y = 0 or x =—E/a corresponding to the classical turning point as shown in Figure 
Dele 


y=- (2.38) 




















urning point__|LE 

















x 


U(x) =-ax 


Figure 2.7: _ Classical turning point atx =—E/aand y=E. 


Using the solutions of the Airy equation given by (2.22) and replacing the 


variable y , the solution in terms of the original coordinate x, the energy E, the con- 


stants a and fi, the particle mass m and the constants C and D can be written as: 


E 
a 
~ fe Se =. 
h h 
) (*) 


The constants C and D are determined by the boundary conditions and for 


W(x) = CAi +DBi (2.39) 




















the above potential D = 0 since Biis exponentially growing in the barrier region or nega- 
tive x direction. 
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c. Square Quantum Well Under an Electric Field 


Generally, when an electric field is applied to a quantum well structure as 


schematically illustrated in Figure 2.8, the profile of the potential will be changed [1, 12]. 
V(x) 





Figure 2.8: | Single quantum well structure with forward applied bias. 


This change is given by the equation: 
V(x, F)=V(x,0)-exF, (2.40) 
where V(x,0) is the potential profile of the quantum well, F is the applied electric field 
in V/m,e is the electron charge and x is the associated spatial coordinate. 


By substitution of Equation (2.40) into the Schrédinger equation we arrive 


at the following formula: 





iv OF) 
in, be +(V,-exF )(x)¥ ,(x) = EY ,(x). (2.41) 


. . . . . “th 
Since the potential energy V,is a constant within the j" layer of the quan- 


tum well, the solution to Equation (2.41) can be written as a linear combination of Airy 


functions [16] as: 


¥(x) = C,Ai(p,(x))+ D,Bi(p,(x)), (2.42) 
and p,(x) is defined as 


16 








—-Fx-7, 
Pi) =7—Gas> (2.43) 
iv F? 
ES 
where m’ is the effective mass associated with the j’" layer and 7 ; 18 given by: 
E-V. 
nj -| : } (2.44) 
e 


where F is the energy of the electron. 


The boundary conditions given by Equations (2.7) and (2.8) are still appli- 
cable to the Airy functions solution, given by Equation (2.42). An arbitrary coordinate 
system was chosen to be at the upper left end of the potential well, for simplicity pur- 


poses, as we illustrated in Figure 2.8. 


The boundary conditions between two neighbor layers jand j+lof the 


quantum well, give the following relations between the coefficients: 


C; = 

—_ (2.45) 
Cia {Ai(B)..)Bi (B))- Ai (B,..) BUB,)} + Dy. {Bi(B,..) Bi (B,) - BiB, )Bi (B,,..)} 

D; a 

n (2.46) 


Cys {E AUB AI (Bi) — Al (BAB) + Dis {EAB Bi (B,,.)- Ai (B, BIB...) 


where ,,,and are defined as: 


jtl 








a (2.47) 
2m ef 
B, = ( x; 7;) Kr ’ 
and ¢, is the ratio between the effective masses of the two layers. 
. 2/3 

mM, 
o,= ea ; (2.48) 

sig +1 
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Equations (2.45) and (2.46) can be written in matrix form as: 


oy 
D; (2.49) 


7 Ai(B,,,)Bi (B,)—€,Ai (B)..) BB) Bi(B,,,)Bi (B;)— ¢,Bi(B, )Bi (B).) \[ Cin 
,Ai(B)Ai(B).)-Ai(B)AUB,..) FAB) Bi (B,..)— Ai (B,) BiB.) )\ Dia) 


A repeated application of boundary conditions, similar to that of the unbi- 


ased case, the coefficients of the outermost layers can be related as: 


C C C 
|= oa Ml, | "(mm ‘| “| (2.50) 
D, Dy My, My }| Dy 


In the case of bound states, the wavefunctions in the outermost layers 


should decay, which implies that the coefficients that make them grow (D, andC,, ) 


should equate to zero. Thus, Equation (2.50) reduces to: 


0 M>, My, Dy 
This implies that, 
m,,(E) =0, (2.52) 


which corresponds to quantized energy states E inside the well under the electric field. 


Figure 2.9 represents a single Al, ,Ga,,As quantum well of 40 Angstroms 


width where the applied field is2x10’ V/m. We represent the potential of the well using 


a blue line, the wavefunction with the black line, and the energy level with the red line. 
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Figure 2.9: | Single quantum well structure with large forward applied bias 
(2x10' V/m). 


For example, the bound state energies in the quantum well were obtained 


by plotting m,,(£) as a function of energy, as shown in Figure 2.10, and identifying the 
energies where m,,(E) =0. For the parameters used, there was one bound state in the 


quantum well, at energy of about0.089 eV. The extent of the wavefunction beyond the 
right side of the barrier indicates tunneling through the barrier due its finite width. This 


is responsible for leakage current in quantum well detectors at low temperatures. 
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Graph of matrix element a,, Vs Energy E 
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Figure 2.10: Zeros of the matrix element m,, determines the energy characteristic 
values. 


E. SUMMARY 


In this chapter, we developed the necessary tools to comprehend the solution of 
the Schrédinger equation in a quantum well structure, with and without the presence of 
electric field. The use of transfer matrix method with Airy functions enabled the analysis 
of quantum well structures under external bias. In the next chapter, response of a QWIP 


made from Jn,Ga,,As/GaAs step quantum well will be simulated and compared with 


available experimental data. Ideas such as continuous states and “normalization” of them, 


Fermi’s golden rule for optical transitions and probability current will be developed. 
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Hl. IMULATION ANALYSIS AND DATA 


A. INTRODUCTION 


The simulation program developed in Chapter II was employed to analyze a 
QWIP detector designed by Kevin Lantz [25] and experimentally studied by Michael 
Touse [26], and Yeo Hwee Tiong [27]. The comparison between the analysis and ex- 
perimental results is focused on the infrared absorption and detector responsivity. 


B. DESCRIPTION OF QWIP SAMPLE 


The fabricated step QWIP consists of a multiple quantum well structure. The 


layer specification of the complete detector is presented in Table 1. 


Thickness | Concentration 
Substance} Mole % (A) 


GaAs | 5000 1E+18 





SE 
=e ae 
if ee 
et 
ss 


Table 1. Specifications of the test QWIP (From [26].) 


To increase the active layer thickness of the detector, 25 step In,Ga,_,As quantum 


wells were employed. A single period of the step quantum well is schematically shown in 
Figure 3.1.The width of the well, and the step, are 40 Angstroms each, while a 300- 
Angstrom barrier of GaAs was placed between, in order to reduce the tunneling losses. 
The sample was grown by molecular beam epitaxy (MBE), which offers high accuracy 


for the design parameters, such as doping and layer thickness, since we have disposition 


of one or more pure materials onto the crystal wafer, one layer of atoms at the time, under 
ultra-high vacuum. The preliminary design of the QWIP predicted an IR absorption peak 
at 10.2 um [26]. 
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GaAs 





Figure 3.1: | Diagram of the designed and simulated QWIP (From [26].) 


C. DESCRIPTION OF THE COMPUTER MODEL 
The structure used in simulation is a In,Ga,_As single step quantum well as 


shown in Figure 3.2. 
0.25 
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Figure 3.2: Potential diagram of simulated Jn,Ga,_,As quantum well structure. 
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This quantum well is divided into four regions. The percentage of Jn, in the first, 


and the fourth regions is zero, while for regions two and three is 30 and 10 percent, re- 
spectively. Similarly the width of the regions one and four is considered infinite, while 
regions two and three are 40 Angstroms wide each. 


The behavior of a In,Ga,_,As single quantum well was investigated in the three 
main cases. In the first case we have studied the In,Ga,_,As quantum well without the 


application of an electric field, in the second case the electric field is present, and finally 
in the third case the electric field is present but its value is relatively small, which needed 
a special treatment. The main difference between these three cases, is that, in the first 
case, we have used the exponential solution of the Schrédinger equation, while in the 
case of large and small bias we took advantage of the Airy and asymptotic Airy solutions, 
respectively (Appendix C). 

1. Simulation Parameters 


For our numerical calculations, the values of the effective mass m’ , potential 


height V , and energy gap E, for the In,Ga,_,As quantum well, were employed as: 





eal (3.1) 
m= F; 7 . 
is (0.067m, ) 
0.028m, 


where m, is the rest mass of the electron and ”x ; 18 the mole fraction of In in the region 
j . The bandgap energy E, (ineV ) of the four layers and the corresponding barrier 


heights V, can be obtained by the following formulas according to [1]: 


E} = Ef =1.519+1.247(“x,), (3.2) 
E? = E; =1.519-1.102(%x,), (3.3) 
V, = 0.62(E/ -E*). (3.4) 


Equations (3.1) through (3.3) characterize the physical, structural and electronic 


properties of Jn,Ga,_,As under specific conditions. 
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D. STEP QUANTUM WELL WITHOUT AN APPLIED ELECTRIC FIELD 


if Bound State Energies 
The bound state energies are calculated using the transfer matrix method (TMM), 


according to the derivation of Subsection E-2 of the Chapter II. The energy eigenvalues 


are derived by the condition m,, (E) =0, given in Equation (2.14). In order to find the 
values of E that satisfy the above condition, m,, (£) was evaluated for energies from the 
bottom of the quantum well (EZ =0), to the top (E =V, =V,), using small increments AE . 
In the numerical analysis, the m,, (E ) derived at each energy was multiplied by the corre- 


sponding value at the previous energy and, if the resulting quantity is negative, then 


Mm, (E£) had crossed a zero. This procedure allowed us to determine the bound state en- 
ergies in the quantum well. 


In the case of the step quantum well, there are four layers, which require three 


matrices to relate the coefficients of the outermost layers as given in Equations (3.5) to 


(3.7): 
C 6: 0 m, m., \(0 m.,D 
Jem a}em(s)e(m a )-{ WD a 
2 1 1 M,, My 1 m,,D, 
= iE 5 = 1h 4 IVE | = = ) . 
D, D, D,) \my mz )\D,) \m3D 
C, CO; C, m,, | Mm," 0 m,,'D, 
(F Jeae Jeanaeae (3 }=[ ait a D =| 301 ; (3.7) 
4 3 1 M), M,, 1 mM,, D, 


1 21 321 
uv? May ? May ? 


v 





where M,,j=1,2,3, is given by Equation (2.12) and m where u,v =1,2, 


are the matrix elements of M, ,the resulting elements of the multiplication of M,M,, and 


the resulting elements of the multiplication of M,M,M,, respectively. 


The next step involves normalization of the bound wavefunction, which allows 


the determination of D,. Normalization implies: 


[Yor dx=1, (3.8) 
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where ‘¥(x) and ¥"(x) are the wavefunction and its complex conjugate, respectively. 


For the parameters used for the step quantum well, there was only one bound state 


in the well as shown in Figure 3.3, alone with the normalized wave function. 


0.25 








Potential Height (eV) 














X (Angstroms) 


Figure 3.3: In.Ga,_,As step quantum well ground state energy without applied bias. 


2: Continuum State Energies 

The responsivity of QWIP detectors is strongly dependent on the location of the 
excited state [1, 17, 18]. To avoid the suppression of photocurrent due to tunneling, the 
excited state is usually aligned with the barrier which also minimizes the thermionic 
emission of ground state electrons [1]. It is possible to calculate the continuum state 
wavefunctions by using the transfer matrix method, with proper boundary conditions. If 
the electron is incident from the left side of the potential well, as shown in Figure 3.2, 


then the coefficient D, = 0, since it is not possible to have a left moving wave in the right 


side of the well. 


Using the transfer matrix method, the coefficients of the fourth and the first re- 


gions can be related as: 
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C C 321 321 C 
*)=aat| (me | i] (3.9) 
0 D, m, m3; )\D, 
The coefficients D, and C,can be obtained in terms of the matrix elements de- 


fined earlier as: 


321 


m 
et ht 
Dea OS (3.10) 
22 
321,321 321,321 321,321 
_ 321 3217) _ 321 mM, _ mM, Mm, C,—m, Mm; 
C, =m, C,.+m, D, =m, CS ee (3.11) 
22 22 


Similarly, the transfer matrix method gives the coefficients of the continuum 


wavefunctions for the second and third regions: 








mm 
1 
C C 1 1 C, mC, — al 
a>) kel ee sat) | = a (3.12) 
— 1 = M1 — ° . 
D D jl. il _ My 1321 
2 1 21 22 321 WhGos MyM), 
My 211 321 
22 
21 321 
m 
21 12!) 
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Figure 3.4 shows the dependence of coefficients in different layers as a function 


of energy, for the step quantum well structure. 
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Figure 3.4: Coefficients of the continuous states vs. energy difference. 


Normalization of the continuum wavefunctions cannot be executed using Equa- 
tion (3.8), since these wavefunctions extend to infinity. The normalization of the contin- 
uum wavefunctions are usually done either using a hypothetical box around the quantum 
well [18], or using momentum space 6 function normalization as described in [20]. Fig- 
ure 3.5 illustrates 100 continuum wavefunctions up to 0.285 eV above the ground state 
using an arbitrary scale. The red curve in Figure 3.5, illustrates the ground state wave- 
function. The oscillatory behavior of continuum wavefunctions represents the free mo- 


tion of in the continuum. 
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Figure 3.5: | Continuum wavefunctions (blue) and ground state wavefunction (red). 


E. STEP QUANTUM WELL STRUCTURE UNDER AN APPLIED ELEC- 
TRIC FIELD 


1. General 

As described in Chapter II, it is necessary to use Airy functions to represent wave- 
functions when the quantum well is under bias. The implementation of Airy’s function 
routines in a computer proved to be a relatively complicated task, especially when the 
electric field is small. Thus, a “weak field limit” was placed where the built-in Matlab 
Airy functions give numerical overflows, so lower biases than this threshold were treated 
using the asymptotic expressions of Airy functions. This limit was found to be about 
10°V/m. The analytic coding used for this case, along with the previously examined 
cases is presented in the Appendix A. 


2 Bound State Energies 


The following shows the calculation of bound state energies in the step quantum 


well under a moderate electric field of 10°V/m. The energy eigenvalues can be found by 
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satisfying the condition (2.52), while the coefficients of Airy functions for different re- 
gions can be obtained as: 
Cin, (.=7..6 =m. =0; (3.14) 


D,=0,D, =m, ,.D,=n0..D, =1; (3.15) 


3 23 123 = : 
where m’,,,m,,,,M,,, #,V =1,2, are elements of M,,M,M,,M,M,M, matrices, respec- 


tively . There was only one bound state in the well as shown in Figure 3.6: 
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Figure 3.6: | Ground energy and wavefunction in the step quantum well under an elec- 
tric field of 10°V/m. 


In the case of negative applied electric field of -10° V/m the ground state energy 


increases compared to the unbiased quantum well as shown in Figure 3.7. 
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Figure 3.7: | Ground energy and wavefunction in a /n,Ga,_,As quantum well under 
—10°V/m applied field. 


It is well known that the bound state energies in a quantum well do not shift when 
measured from the center of the well [1, 17]. Since our origin is at the edge of the well, 


the amount of the energy shift can be estimated using: 


le| FL 
shift =— ot (3.16) 





where eis the electron charge, F' is the applied bias, L is the length of the well which 
confines the energy state. This implies that for a 10° V/m applied field, the shift is about 
—0.002 eV from the unbiased position, and for —10° V/m the shift is about 0.002 eV. 
This gives an energy separation of 0.004 eV for +10° V/m field which is very close to the 
simulated value of 0.005 eV. 

In the case of application of an electric field of 7x10* V/m, the asymptotic 


expressions of Airy functions were needed to obtain the ground state energy and the 


associated wavefunction are illustrated in Figure 3.8. 
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Ground energy and wavefunction in a quantum well under7x10* V/m. 


Figure 3.8: 


In the case of the same bias applied in the negative direction, the program gives 


that the ground energy E, lies at0.091 eV, as shown in Figure 3.9. 
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Ground energy and wavefunction in a quantum well under—7x10*V/m. 


Figure 3.9: 


a1 


Again, the difference between the values of the ground energy level can be exp- 
lained by Equation (3.18) since, for an 7x10* V/m applied field, the shift is about 
~7x10~ eV from the unbiased position, and for -7x10* V/m is about7x10~ eV. This 
gives an energy separation of 0.00143 eV for +7x10* V/m field which is very close to 
the simulated value of 0.00137 eV. Finally the lowest bias that we have achieved in our 


simulation was +7 x10* V/m or+0.7 kV/cm, which is much smaller than typically used 
in QWIP detector operation (10 kV/cm) [1, 16]. 


3. Continuum State Energies 


The continuum state wavefunctions are obtained by assuming the electron is inci- 
dent from the right side and setting the following boundary conditions for the coefficients 
of the outermost layers: 

Di 0 (3.17) 

Then applying the TMM, the coefficients in each layer can be obtained as: 


123 
My, 
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C,=m,C,+m,;D,=m,,C,+m,,; L-( a , (3.21) 
22 22 
123 123. 123 123.123 
m mm, C,—-m5m 
a 193 1237, _ ...123 123 My, _ | My My © 12 My) 
C, =m, C, +m, D, =m C,-—m. ay = 123 ‘ (3.22) 
22 My, 


For the normalization of the continuous wavefunctions, the technique described in 
the Appendix of Ref. [21] was used, since 6 function normalization was not valid due to 


spatial dependence of the potential energy. 
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The absolute values of the coefficients as a function of continuum energy are 
illustrated in Figure 3.10. In the same figure, we observe the oscillatory behavior of the 
coefficients due to the resonances formed in the continuum region, due to constructive 


interference of incident and reflected waves. 
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Figure 3.10: Absolute value of coefficients under 10° V/m forward bias. 


The first 100 continuum wavefunctions of the tilted well are illustrated in Figure 


3.11, along with the ground state wavefunction. 
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Figure 3.11: Continuum wavefunctions (blue) and ground state wavefunction (red) for 
the step quantum well under 10° V/m forward bias. 


F. BOUND-TO-CONTINUUM TRANSITIONS IN A QUANTUM WELL 


1. Fermi’s Golden Rule 
It was shown by Enrico Fermi [9], that the transition probability W of a quantum 
system from an initial state / to a final state F is given by: 


2 


Woe = = tz VV, |,)| 6(£, -E, -ha), (3.23) 
F 





") 





where Y,, and Y, are the final and the initial wavefunctions, respectively, FE, and E, 
the associated energies, V,, is the interaction potential and fw the interaction photon 


energy. The photon interaction potential in dipole approximation is given by [2]: 


1/2 
Th 
Vs = “(4 EP. ’ (3.24) 


m \ 26 €@c 
where ¢ is the electron charge, m’ the electron effective mass, I , is the incident photon 


flux, ¢,andé characterize the electric permittivity, c is the speed of light, @ is the 
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frequency of the photon, p, defines the momentum of the electron and € is the 


polarization vector of the photon. 


Ze Oscillator Strength 


The oscillator strength f,,. for a transition from an initial state to a final state is 


defined as [1]: 
2m 2m 2 
Tee =<" (E; ~E,)(zy == (Er -E,) |. 7 (z) oF ,(2z) dz, (3.25) 


where E,, and Y’,,are the final energy and wavefunction, respectively, E, and Y, are 


the initial energy and wavefunction, m” is the effective mass, / is Plank’s constant 
and z is the direction of growth. 
Generally the oscillator strength is an indication of how strong the transition is, 
which directly translates to the performance of a QWIP detector. 
Figure 3.12 shows the calculated oscillator strength for the step quantum well, 
under zero bias, using the wavefunctions derived in Subsection D-2. The energy was 


measured from the ground state. 
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Figure 3.12: Oscillator strength versus energy for bound-to-continuum transitions in an 
unbiased In,Ga,_,As step quantum well. 
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Under an applied electric field the oscillator strength can be calculated using the 


wavefunctions derived in Subsection D-3. Figure 3.13 shows the oscillator strength for 


10° V/m applied bias. 
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Figure 3.13: Oscillator strength under 10° V/mapplied bias vs. energy. 


The oscillatory behavior of the oscillator strength is due to the constructive inter- 
ference of incident and reflected electron waves, from the tilted left quantum well barrier. 

This oscillatory behavior becomes denser, as the strength of the electric field was 
reduced, as illustrated in Figure 3.14. As the bias approaches zero, the oscillations be- 


come closer and closer, to form a continuous curve as shown in Figure 3.12. 
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Figure 3.14: Oscillator strength under 7x10* V/m applied bias vs. energy. 


The origin of oscillations as a function of electric field can be understood by esti- 


mating the amount of energy change (AE ) required for obtaining path difference between 
incident and reflected electron waves (Az) by one wavelength. This can be estimated by 


using the linear potential generated by the electric field as: 


AE =|e| FAz. (3.26) 


It can be seen from Equation (3.26) that as the electric field becomes smaller, the 
amount of energy change needed to get one wavelength path difference reduces. This 
implies that the constructive interference occurs at a rapid rate, as the energy of the elec- 
tron goes up, resulting a rapid oscillation of the oscillator strength. 


3. Absorption Coefficient 


The performance of photodetectors can be conveniently characterized using the 


absorption coefficient which is defined as: 


a7 


ho x (number of transitions per unit volume and time) 





a(ho) = (3.34) 


incident energy flux 


For a QWIP detector the absorption coefficient can be expressed using the density 


(, 


where JN, is the doping density in the well, Lis the total length of the doped region, m, 


of states and oscillator strength as [1, 2]: 


a(ha) = 





a 
dx 





(3.35) 








Ne2h?L | 2m; 
(nw) (E, -V,) 








2(m") ee 


and V, are the barrier effective mass and potential height, respectively, and fa is the 


incident photon energy. Figure 3.15 shows the oscillator strength and the density of 


states for the step quantum well under a 10° V/m electric field. 
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Figure 3.15: Oscillator strength and density of states vs. energy. 


The absorption coefficient was calculated for the step quantum well assuming the 


doping density to be 10'* cm“ and the result is shown in Figure 3.16. The calculated 
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absorption coefficient is in good agreement with that typically observed experimentally 


(800 cm“) [22]. 














1400 





1200} ------i------ 


{ 
— 
fo} 
fo} 
oO 


1) ee eee neeeerepee sy cunereeeeeee Spo eee 


G00 saceudsacads 


Absorption coefficient (cm 1) 


400} - -----i------ 


7.5 Rees ERE Sere Serene oc meres oreemeree: Game? veen Scere serS 4 











Photon Energy (eV) 


Figure 3.16: Absorption coefficient vs. photon energy under 10° V/m applied bias. 
G. COMPARISON WITH EXPERIMENT 


The step quantum well structure used for the simulation in the present Chapter 
was fabricated and characterized by Touse [26] and Yeo [27] using photocurrent spec- 
troscopy. In this section, the measured responsivity of the detector is compared with the 
simulated responsivity using the absorption coefficient and outgoing probability current 
of the excited electron. 

1. Responsivity 

For a detector the normalized voltage or current signal per incident optical power 


is called responsivity [24]. We can express the responsivity R as: 
R=n—g, (3.36) 


where 77 is the quantum efficiency, e is the electron charge, hv is the average photon 


energy and g is the photoconductive gain. 
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Taking into account all the geometrical and optical factors, the test detector re- 


sponsivity can be estimated using [23]: 


2-I,,,-R_,-d-Sw 
pw tsa 2 toa Regd Sw (3.37) 
Pret V, * Aut Ze T. 


ref ZnSe GaAs 


where /,,, refers to the photocurrent, P,., is the optical power incident on the test detector 


det 


power, R 


ee 1S the responsivity of the reference detector, d- Sw refers to the reference de- 


tector surface area, V 


ve 1S the generated by the reference detector voltage, A,.,is our own 


and 7, , are the transmittances of Zinc Selenide and Gallium 


detector surface area, T, oe 


ZnSe 


Arsenide, respectively, of the cold head where the detector was situated. 


During the experiment the temperature was kept fixed at 40K, so the photocur- 
rents for the test detector under different applied electric fields were measured, and the 


responsivities calculated using Equation (3.37). 


The experimental responsivity as a function of wavelength, for the step quantum 
well detector, for a set of biases is illustrated in Figure 3.17. The shifting of the peak po- 
sition in Figure 3.17 is mainly due to the linear Stark effect in step quantum well as ex- 


plained in Chapter II. 
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Figure 3.17: Experimental responsivity vs. wavelength atT = 40 K (From [27].) 
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2. Absorption 

The normalized absorption coefficients calculated using the approach described in 
Subsection F-3 of the present Chapter, as a function of wavelength for the bias voltages 
used in the experiment are shown in Figure 3.18. The shifting of the peak position can be 


clearly seen from Figure 3.18. 
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Figure 3.18: Computer simulation of absorption vs. wavelength for the step quantum 
well structure. 
3. Photocurrent 
The photocurrent density can be estimated by taking the product of quantum effi- 
ciency with the outgoing current density of the excited electron. The current density as- 
sociated with a quantum state can be calculated using the wavefunction [24], as: 
hie «ad? yy 
aah So cle ad : (3.38) 
2mi dx dx 
where fi is the reduced Plank’s constant, m is the electron mass, x 1s the direction 


where the current is evaluated and V is the wavefunction. 


Assuming the wavefunction has the following form: 
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CAi +DBi, (3.39) 


and after some algebraic effort, the quantity in the bracket of Equation (3.38) can be ob- 


tained as: 


wy" dv iy d¥ 
dx dx 





=(C’D- D°C)( AiBi — Ai Bi) (3.40) 


Using the Wronskian (W) of the two Airy functions [40] the second bracket in 
Equation (3.40) can be further simplified as: 
Ai Bi 


he nd By eee ree (3.41) 
Al Bi 


5 








The current density given by Equation (3.38) can be expressed using the result in 
Equation (3.41) as: 
h * * 
eal G D-D'Cc), (3.43) 
271 m| B | 
where £ is given by Equation (2.20), and appears in the denominator of Equation (3.43) 


due to the argument of the Airy functions of Equation (3.39),which according to Equation 


(2.39) is: 
tars 
x+— 
a 


Al 


(3.44) 


Equation (3.43) represents the total current density which should be zero, since 
the probability density is time independent [20]. This requires that the coefficients 


C and D must be related as: 
C*xD =D*xC. (3.45) 


In our simulations, it was found that this condition holds for all the regions of the 
quantum well structure. This indicates that the current density given in Equation (3.43) 
consists of two identical components, (a) outgoing current and (b) incoming current, 


which cancel to give zero net current. 
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In order to estimate the photocurrent after excitation of the electron from the 
ground state, it is necessary to extract the outgoing component from the current density 
that is given in Equation (3.43). The incoming and outgoing components of the current 
are not obvious in the case of the linear combinations of Airy functions, while they are 


usually obvious for exponential wavefunctions obtained for the zero bias case. 


However, it can be shown that the linear combinations shown in Equations (3.46) 
and (3.47) produce pure outgoing and incoming currents by substituting in Equation 
(3.38): 

Y , = Ai-iBi, (3.46) 


Y= Ai+iBi. (3.47) 


The magnitudes of current density corresponding to the outgoing and incoming 


wavefunctions are given by: 





Jug = av ai 3.48) 
2mi dx 

i We, —WY os : (3.49) 
2mi dx 


Under a positive bias, the photoexcited electrons moves to the right of the quan- 
tum well and the wavefunction in outermost right region, given in Equation (3.50) can be 
used to determine the outgoing current. 

Y,=C,Ai+D,Bi, (3.50) 
where C,and D, are the coefficients defined in Subsection D-2 of Chapter I . 

Equation (3.50) can be rewritten using the outgoing and incoming wavefunctions 
given in Equations (3.46) and (3.47) as: 

WHO ¥ 4D Px, (3.51) 


where the constants C and D are related to the coefficients C,and D, ,since: 


cles (3.52) 
2 
iD 
p=, (3.53) 
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The corresponding currents are given by: 


2 























See ear aa (3.54) 
am| B | 2 
‘ 2 
| ba (3.55) 
am| B | Z 
The photocurrent density was calculated using the outgoing current given in 
Equation (3.54) and the quantum efficiency 7 = a@Las: 
: 2 
jag (3.56) 
mm | B | 2 








where @ is the absorption coefficient and L is the total length of the quantum well struc- 


ture. 


Figure 3.19 shows the calculated photocurrent for the biases used in the experi- 


mental measurement of the step quantum well detector. 
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Figure 3.19: | Calculated photocurrent vs. wavelength. 
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Figure 3.20 compares the absorption and normalized photocurrent density with 
the measured normalized responsivity at 0.81 V forward bias. It can be observed from 
this figure that the simulated photocurrent density provides a good description of the ex- 
perimental observations. The absorption is found to be higher than the photocurrent den- 
sity in the longer wavelength regime. This is mainly due to the difficulty of extraction of 
electrons after excitation due to the triangular barrier on the right side of the biased quan- 


tum well. 
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Figure 3.20: Measured responsivity, simulation photocurrent density and absorption vs. 
photon wavelength for 0.81 V forward bias. 
H. SUMMARY 


In this chapter, we presented the calculation of absorption and photocurrent den- 
sity for the step quantum well detector under an external bias. The calculated results 
were compared with the experimental data and found to provide a good agreement which 
validated the accuracy of the model employed. In the following chapter, we will present 
a simulation code based in EKF for angle only tracking using an IR sensor. Ideas such as 
initialization using two synchronous bearing and EKF angle tracking will be developed 


and finally, selected results will be presented. 
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IV. PERFORMANCE OF QUANTUM WELL INFRARED 
PHOTODETECTORS IN TARGET TRACKING 


A. INTRODUCTION 

Quantum well infrared photodetectors constitute a tool for target tracking since 
they offer the advantage of passivity. As we have discussed in the previous chapters, tak- 
ing advantage of our ability to design multi-wavelength tunable infrared detectors, we 
can combine them with batch processing or recursive type algorithms in order to track 


targets passively. 


In this chapter we simulate in Matlab the scenario of bearing-only tracking of a 
target using a recursive algorithm based on the Kalman filter. We investigate the bearing- 
only tracking problem in two main cases, namely tracking a non-maneuvering target and 
tracking a maneuvering target using either one or two sensors. Since we study only the 
tracking problem, we suppose that the detectors are tuned in such way that they provide 
the highest contrast and avoid the extensive clutter so they generate relatively reliable 
bearings. More analytic information about the coding used can be found in Appendix B. 
Finally the BOT problem is considered trivial since limited information is available 


B. DESCRIPTION OF THE MODEL 


The simulation model and the encounter geometry are illustrated in Figure 4.1. 


» 





Sensor 1(0,0) Sensor 2 (10km,0) x 


Figure 4.1: | Encounter geometry of the Sensors and the target. 
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Sensor | is located at position (0,0) of the local coordinate system, while sensor 


2 was placed at (1 Okm, 0) . These sensors, since they are IR, can make measurements of 


the bearing of the target only. The first step of the bearing-only tracking using an ex- 
tended Kalman filter is the initialization of the filter. Initialization in the case of EKF is 
sufficient if we estimate two positions of the target. This is possible either by taking ad- 
vantage of the observer and target relative trajectory in the case that we use one sensor, or 
by assuming synchronous multiplatform bearings. BOT using asynchronous bearings will 
not be examined. So without the loss of generality in our case, we assume two stationary 
IR observers and synchronous initialization and tracking bearings. 


t. INITIAL ESTIMATE OF THE TARGET POSITION 


The initial estimate of the target position is determined using two sets of time syn- 
chronous bearings. We have used the fact that, in order to initialize our filter, we do not 
need to have the exact position of the target, only a reliable estimate of its position. The 


simplified geometrical problem is illustrated in Figure 4.2. 






ne, 
— 

BOX. Yds. 
t 


Ul 


O Sensor 1(0,0) A Sensor 2 (10km,0) 


Figure 4.2: | Encounter geometry of the initialization problem. 


We suppose that the first set of the synchronous bearings from our IR sensors are 


6, and @, as shown in Figure 4.2. The reliability of the bearings depends on the quality of 
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the sensors, the sensors-target geometry and the distance between the sensors and the tar- 
get. Since we examine the general case of the problem of initialization, we suppose that 
the sensors we use give a random error in bearing of maximum magnitude of approxi- 
mately 2°. This error is considered very high for sensors designed to provide only the 


bearing, and it can be positive or negative. 


So, examining sensor 1, the angular real position of the target must be between: 





OS SOR (4.1) 

where 
6- = 6, —|rand(0,1)-29, (4.2) 
6" = 8, +|rand(0,1)-2°|, (4.3) 





while bearings 0 and @ are shown in Figure 4.2 with the blue dashed line. 


Treating the sensor 2 case, using the methodology used for sensor 1, we observe 
that the target must be inside the quadrilateral BCDK . In order to find the coordinates of 
BCDK we apply Euclidian geometry to each triangle of Figure 4.2. So, examining trian- 
gle OBA and applying the sine law, we get: 

sin(@*) _ sin(z—@;) _ sin(@) 
AB OB OA 





; (4.4) 


where AB, OB, OA are the lengths of the sides of the triangle as shown in Figure 4.2 
while @} refers to the third angle of the triangle, namely angle OBA. From Equation 


(4.4) we can find the length OB as: 





C5204 = (4.5) 
sin(@}) 
where OA is the distance of the two sensors while angle @} is obtained since: 
6} =2-0* —-(a7-67) = 6; -O. (4.6) 
Obtaining the length of OB we can find the coordinates using: 
PA Tee Ch eer ab a (4.7) 


sin(O; — 05) 
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sin(@; )sin(@;*) 


= OBsin(60*)=OA 
vi ( ‘ie sin(O; —0*) 


(4.8) 

We observe that finding the coordinates of point B involves only the sensor dis- 
tance and the uncertainty bearings obtained by the two sensors. Following the same pro- 
cedure the coordinates of point C are given by: 


sin(@;) cos(@*) 








=0OCcos(0*)=OA ; 4.9 

x, (O") sin(@; —67) (4.9) 

“Sonsini =on  ). (4.10) 
sin(@; — 6") 

For the coordinates of point D our calculations give: 

~0neeoj=04 (4.11) 
sin(@; —6,) 

22 OD j 20s (4.12) 


sin(@; —6,) 
while for point K 


sin(@; )cos(@, ) 


=OK 0-)=OA 
Le cos(@ ) sin(O; -0-) 


’ (4.13) 





sin(@; )sin(@, ) 


= OK sin(@,)=OA 
ue a sin(O; —0-) 


(4.14) 


The next step involves determination of the maximum between the two diagonal 
lengths BD andCK. This determination is used for initialization of the covariance ma- 


trix P of the extended Kalman filter. The determination is illustrated in Figure 4.3. 


(0, 4) 





(x,,0) 


Figure 4.3: | Encounter geometry for the determination of the diagonal lengths. 
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The diagonal length can be extracted either using the angle @ of the orthogonal 


triangle MCK whose inverse tangent is: 


g=tan" pas — ; (4.15) 
Yo V4 
while the length CK according to this approach: 
X,—-x. 
CK = l=) (4.16) 
sin| tan"! (x8) 
Yo Na 
or using the well known approach for finding length: 
BD= (x,-x,) +(¥,-y,) ; (4.17) 


Since the target can be anywhere inside the quadrilateral we assume without the 
loss of generality that it is located at: 


ete ty +% 





, 4.18 
; 4 (4.18) 
a cn (4.19) 


We have implemented the above considerations in Matlab and we have produced 


Figure 4.4, for the target initialization: 
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Figure 4.4: Matlab initialization using two sets of synchronous bearings. 
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The position of the target during the first and the second time that the bearings 


took place is illustrated in Figure 4.5: 
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Figure 4.5: The position of the target is located inside the uncertainty bearing box. 


D. EXTENDED KALMAN FILTERING CONSIDERATIONS 


1. General 
The problem considered in the bearing-only tracking is that the target motion is 
modeled in linear Cartesian coordinates, while the measurements are in polar coordinates. 
We begin by assuming the following non-linear discrete-time system: 
x(k +1) = f (k,x(k))+(k,x(k))y,, 


2(k)=h(k,x(k))+@, ae 


where x provides the target position data, z is the measurement, f andh are non-linear 
vector-valued functions, g is a non-linear matrix valued function and v,, @, are uncor- 
related Gaussian processes. 
According to [28-30], it can be proved that the prediction using an extended Kal- 
man filter is given by: 
Resiv =F (ete) (4.21) 
Prosi = FP yl +G,Q,G;, (4.22) 
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where X,,,, is the state vector at time k +1 given its value at time k, F,,,, is the cova- 
riance prediction matrix at time k +1 given its value at time k, F, and F; are the gradient 
and the transpose gradient of f (k, Kup )s G, and G/ are the zero order term for g (k, ae 
evaluated at <,,, and its transpose, respectively and, finally, Q, is a matrix related to the 
noise inserted in the system. 

We define Q, according to [27] as: 

0, = 70, (4.23) 
where Q, is the sampling interval matrix and q arbitrary chosen parameter of the design 
[30]. 

The update of the measurement according to [28, 29], can be expressed as: 


resi = Kes + Ky, (aa = h(k Fl, Kessit )) > (4.24) 
Pees = (J ~K, Ay ) Poa (J ~— Ky, Ay, ) +K, A Kia» (4.25) 


where z,,, 1s the measurement at time k +1, / stands for the associated to the problem 


identity matrix, H,,, is the gradient of h(k +1, 5.15 ) and K,,, 1s defined as: 
-l 
Kya = Posie Lin ze eee ae 25 Russ ) 2 (4.26) 


where R,,, 1S associated to the uncertainty of the measurement at time k +1. 


2. Bearing Tracking Expansion of EKF 
The initialization in the EKF BOT case provides the vector of the two positions of 


the target after the second set of the synchronous bearings. This vector is: 
T 
x= (x, pes y,) ; (4.27) 
where x,, y,,x,, y, are the estimated possible positions of the target obtained during the 


initialization process by the two sets of the tautochronous bearings. Since our state vector 


is in the form: 


f=(x v, y V, Ie (4.28) 
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where x, y is the position of the target and v, ,v, the associated velocity, we use matrix D 


to obtain: 
xX — Dx 5 (4.29) 
where the matrix D is defined as: 
1 0 0 0 
l/dt O- -Il/dt 0 
D= , (4.30) 
0 1 0 0 


0 1/dt 0 —1/dt 


where dt is the time between the measurements. 


The initialization of the covariance matrix P involves the maximum diagonals of 


the quadrilaterals and it is given according to [29, 30] by: 


ad 0 0 0 
p (0 @ 0 0 rer 
initial 0 0 d? 0 9 . 
0 0 0 @& 


where d, is the maximum distance of the diagonals of the quadrilateral formed by the 
first set of bearings, while d, refers to the second. During the EKF update and for the 


rest of the calculation the covariance matrix P is obtained by Equation (4.25). 


The bearing-only measurements are 0 angle measurements, so they can be ex- 


pressed by: 
z(k) =0(x(k)) +a. (4.32) 
where 
O(x(k))= actan(* (4.33) 
Ve 


The gradient of h(x) =6(x) for the bearing-only measurements can be expressed 


by the following matrix: 


4, -E)| Boas Hy ei 0), (4.34) 


Ox e+y? x+y? 








where x and y are obtained by Equation (4.27). 
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3. Maneuvering Target Considerations 


According to [29] and [30], in the maneuvering case target, the equations of 


measurement and update do not change. However, the transition matrices should be cal- 





culated as: 
sin(@A) 1—cos(@A) 
oO o 
oe 0 cos(@A) 0 seu) (4.35) 
0 1—cos(@A) 1 sin(@A) 
oO oO 
0 sin(@A) 0  cos(@A) 


where @ is the turn rate and A is the sampling interval. The determination of F“ is 
possible by application of simple target motion analysis for a rotating mass around a 
complete circle. The turn can be either clockwise or anticlockwise. 
4. Trajectory Error and Measurement-Angle Chi-square Formulation 
The error is determined by the difference between the trajectory that the EKF pre- 
dicts and the real trajectory of the target: 


ee a fs 
Vetter y Yreal 


are the difference between the predicted by the EKF x, y and the real 


(4.36) 


where x 


error ? Nv eior 


coordinates x,..,, ),.,, of the position of the target. As a consequence the total distance 


x 
Total Error = { 
y 


error 
Xerror Yerror r 


error 


real? 


error is obtained by: 





(4.37) 


The Chi-Square values of the angle allow us to observe the angle variations dur- 


ing the measurements. These values according to [27] are obtained by: 


chi? = (z-2) -(H-P-H'+R)"-(z-8), (4.38) 
where z is defined by Equations (4.32) and (4.33), z is the non-linear angular measure- 
ment, H is the gradient defined by Equation (4.34), P is the covariance matrix and R 
involves the error of the measurement. In every sampling interval the chi? value changes 


since the measured angle of the target is different. 
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E. PROBLEM FORMULATION AND RESULTS 


The mathematical considerations above are the basis of tracking a moving target 
both in a straight line and during a turn using either one or two sensors, with or without 
the presence of measurement noise. For the purpose of this thesis we simulated each case 
separately in Matlab, although only selected cases will be presented. Analytic informa- 
tion for the all the cases can be found in Appendix C. In our simulation we assumed that 
the position of the sensors and the bearing measurements were according to Figure 4.2. 
The initialization always takes place using two sensors since these are assumed station- 
ary, but it is possible to use only one if a sensor motion is assumed. In this case the 
mathematical equations should be corrected in such way that the relative trajectory of the 
target-observer is taken under consideration [29]. 


1. Non-Maneuvering Target BOT using One or Two IR Sensors in Noise 
Presence 


In the case of a non-maneuvering target the initialization of the filter takes place 
using both sensors, while during the tracking we can use either one or two without any 


difference in the results. The geometry is illustrated in Figure 4.6. 


x 19° BO Tracking two Sensors, Noise presence,q? = 1000, accel =0g 
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Figure 4.6: True position and track position of straight line moving target using two 
IR sensors in the presence of measurement noise. 
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The error during the tracking and the values of the measurement angle are pre- 


sented in Figure 4.7. 


Track Position Errors: BO Tracking two Sensors, Noise presence,q> = 1000, accel =0g 
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Figure 4.7: Track position error and chi-squared values vs. sample number for non- 
maneuvering target - two IR sensors in the presence of measurement noise. 


In Figure 4.6 we observe that, when we use both sensors for initialization and 
tracking, the filter tracks the flying target almost from the starting point of our simulation 
run until the final point. We have selected a geometry for the simulation which does not 
offer advantage to the EKF since the target flies away from the sensors and the bearing 
angles are not perpendicular. Despite these facts, the filter error generally remains under 
200 m for most of the simulation time, although there is an interval of 40 samples 
where the error increases. 

In Figure 4.7 we also observe that the values of the chi-squared variable remain 
under 9-20 which means that the effort of the filter to correct the bearing angle and, 
hence, to obtain a more accurate position of the target is minimal. This means that the 
filter in unaware that there is a difference between the tracking position and the real posi- 


tion of the target. This result was expected since we are simulating a filter that tracks only 


a7 


with angle measurements, so in a case when the EKF is used with radar sensor that pro- 
vides range information, these effects are eliminated. 
2. BOT of Low Maneuvering Target using IR Sensors in Noise Presence 
In this case the target turn was simulated according to Equation (4.35). We will 
examine two cases in this category. The first case involves a smooth target turn of not 


more than3 g, while in the other case the target will execute a sharp high g turn. The ge- 


ometry for the filter and the target for the first case are shown in Figure 4.8. 
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Figure 4.8: True position and track position of a low g turning target using one IR 
sensor in the presence of measurement noise. 


From Figures 4.7 and 4.8 above we conclude that as long as the target flies on a 
straight course the single sensor filter responds positively. When the target starts its turn 


we observe a slow response time from the EKF which results in the loss of the target for 


almost km .When the filter realizes that a turn has occurred, it turns sharply seeking for 
the true target track. The track subsequently generated by the EKF parallels the true tar- 
get path, matching the true bearing closely but missing in range. So the filter tracks the 
target path only in bearing. This behavior appears because we applied the EKF algorithm 


to track a turning target using only one sensor so, once the sensor loses the real position 
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of a target can only recover the target’s bearing, but not the distance. The error and the 
chi-squared graphs presented in Figure 4.9 supports the previous arguments. 
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Figure 4.9: Track position error and chi-squared values vs. sample number for low g 
maneuvering target - single IR sensor in the presence of measurement noise. 


In Figure 4.9 we observe that using both sensors for initialization and a single 
sensor for tracking results in an increasing error to almost 700 meters before the target 
starts its turn. The filter tracks the target almost from the starting point of our simulation 
run until the final point. The tracking in this case is characterized as low quality since the 
true position of the target after the turn is never obtained. Atsample 28 , as we clearly ob- 
serve, the error starts growing and, at the end of the simulation, it becomes almost 
30 km .We also observe that the values of the chi-squared remain under 0.10 when the 
target is out of the turn loop and close to 0.50 radians when the target is inside the turn. 
The fact that the values are kept low means that the filter assumes that it is tracking the 
target accurately after the turn. This is true in bearing but not in range. Again since we 


combine the EKF with a single angle tracking sensor, the algorithm is unable to regain 
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the real position. The tracking results become worse when the simulation involves higher 


random errors and sharper turns, as we will present. 


To improve the performance in the case of a3 g turn, we employed both sensors in 


the tracking effort. The results are shown in Figures 4.10 and 4.11. 
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Figure 4.10: True position and track position of a low g turning target using two IR 
sensors in the presence of measurement noise. 
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Track Position Errors: BO Tracking two Sensors,Noise presence,q> = 1000, accel =3g 
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Figure 4.11: Track position error and chi-squared values vs. sample number for low g 
maneuvering target — two IR sensors in the presence of measurement noise. 


We observe that using the second sensor has improved the tracking conditions; 
namely, it has reduced the resulting error from 30km to almost 4km. If the simulation 
was free to run for a longer sampling time, the real position and the tracking position 
would eventually converge. At the same time the chi squared values were found almost 
the same. 


3. BOT of High Maneuvering Target using IR Sensors in Noise Pres- 
ence. 


We present two sub-cases in BOT of a high maneuvering target. Firstly, the concept 
of EKF angle tracking using one sensor. In this case, a turn of no more than 9 g is pre- 
sented since, with the particular geometry and with the particular initialization bearings 
and sampling time, higher turns result in a breakdown of the EKF algorithm. So employ- 
ing a 9 g turn of a target, which is moving with velocity 339 m/s, results to the geometry 


of Figure 4.12. 
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x 40BO Tracking single Sensor,Noise presence,q> = 1000, accel =9g 
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Figure 4.12: True position and track position of a high g turning target using one IR 
sensor in the presence of measurement noise. 


We observe that the single sensor combined with the EKF tracks the target relia- 
bly when we have a straight-line motion, but when the turn takes place it follows with a 
time delay and finally manages to recover the correct bearing tracking. The associated 


error and the values of chi-squared are illustrated in Figure 4.13. 
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x 10° BO Tracking single Sensor,Noise presence,q* = 1000, accel =9g 
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Figure 4.13: Track position error and chi-squared values vs. sample number for 
high g maneuvering target — single IR sensor in the presence of measurement noise. 


Figure 4.13 shows that the error is kept small before the target starts the turn, it 


grows to almost 8km during the turn, then it reduces due to the convergence of the real 


and the estimate position but finally it grows due to the divergence of the two tracks. The 
chi-squared values are also higher compared to the other cases due to the high g turn, so 
we observe that the values are over | during the turn, while for the half of the simulation 
time, they are close to 0.1. 

Second, in order to obtain more reliable results we employ the second sensor in 
the tracking effort. We expect the performance of our system to be improved similar to 
the results of the low g turning target. In fact, the tracking results turned out to be better 
than expected, meaning that the IR sensors combined with the EKF algorithm track the 
target with small error. The explanation of the reliable tracking is that in this case the tar- 
get follows a reverse course towards our system. This means that the distance to our sen- 
sor system is reduced, so the suggested rule that the tracking distance should not be much 
higher than the sensor distance is satisfied [30]. The results and the encounter geometry 


are shown in Figures 4.14 and 4.15. 
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x 10° BO Tracking two Sensors,Noise presence,q° = 1000, accel =9g 
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Figure 4.14: True position and track position of a high g turning target using two IR 
sensors in the presence of measurement noise. 
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Figure 4.15: | Track position error and chi-squared values vs. sample number for high g 
maneuvering target — single IR sensor in the presence of measurement noise. 
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Figure 4.15 illustrates that the error remains high during the turn, but it is also 
significantly reduced compared to all the previous cases. In addition, the filter with the 
two sensors recovers quickly after the turn and seeks the real path of the target. The chi- 
squared values are much higher than all the previous cases but they last for shorter sam- 
ple periods, which means that the tracking effort is large during the turn and minimal at 
all the other times. 


F. SUMMARY 


We have investigated a classical algorithm applied to two modern sensors. It is 
obvious that the results are significantly improved in the case where two sensors are used 
for tracking. In real applications more than two sensors may be used at the same time. In 
addition, three dimensional cases are examined in the literature [31]. Other techniques not 
explored in the present thesis may also be used. Some alternative methods involve parti- 


cle filters [29], and combined IR and radar sensors with EKF application [32]. 
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V. CONCLUSIONS AND FUTURE WORK. 


We have modeled the behavior of QWIPs in various conditions of operation 
through a computer algorithm and we have simulated the behavior of an infrared sensor 
in the bearing-only tracking application. Our goal, optimization of the performance of 
QWIPs and their application to target tracking, was achieved by implementing in a com- 
puter fundamental quantum mechanical equations and nonlinear control theory, respec- 
tively. Our conclusions are summarized per chapter along with suggestions for future 
work. 

A. CONCLUSIONS 

In Chapter II we began with the fundamental quantum mechanical concept, the 
Schrédinger equation. By application of this equation in the quantum structures of inter- 
est we studied their behavior. Then, using a well known mathematical technique, the 
method of the transfer matrix (TMM), we were able, not only to have our equations in a 
consistent form, but also to implement them inside the simulation program and study 
them. During the program simulations we were able to observe theoretically predicted 
and experimentally observed concepts, such as tunneling through the quantum well bar- 
rier and leakage current. Finally we developed a set of programming functions using as- 
ymptotic expressions of Airy’s functions to study the behavior of our structure in the case 


of very low biases. 


Chapter III provided a discussion of the fabricated QWIP. A single quantum well 
of the multi-quantum structure was examined in various cases. Namely, we studied the 
behavior of bound, quasibound and continuous states in the case of forward and reverse 
bias, both above and below the arbitrarily placed threshold and in the case of absence of 
electric field. We examined the concept of bound and continuous coefficients of the 
wavefunctions used in the TMM and we developed their value in every region of the well 
and for the bias used. In addition, we examined the performance of the detector by study- 
ing its oscillator strength using characteristic biases and we found the theoretically pre- 
dicted rapid oscillations when we apply very low biases. Moreover, we compared the ex- 
perimental and the simulation results and we concluded that they match. Finally, we de- 


veloped and simulated the concept of photocurrent in the quantum well. 
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Finally, Chapter IV presented the basic concepts of nonlinear extended Kalman 
filtering. We developed the mathematical model of initialization of the filter using two 
sets of synchronous bearings. The algorithm was initialized by the user and, based in that 
geometry, various cases were studied, either using one, or two sensors with and without 
the presence of ambient noise. Results and conclusions were made for the optimal design 
of a bearing only tracking system based on IR sensors. 

B. FUTURE WORK 

There are many concepts based on this work that further research is recom- 
mended. For instance, using the simulation program developed, we can design the desired 
QWIP, predicting its behavior in the laboratory. In addition, we can predict the necessary 
bias to observe the theoretically predicted constructive interference oscillations. Further- 
more, we can design the optimal quantum well to minimize the quantum tunneling and to 


maximize the photocurrent before testing it in the laboratory. 


Finally, it is recommended to expand the EKF simulation program, using the 
same initialization process, so that the number of the sensors exceeds two. In that way we 
can simulate sensors arrays and the expected tracking results will be better. Furthermore 
we can expand the simulation into three dimensional cases. Concluding, it is strongly 
recommended to study similar cases using particle filters and compare the results with the 
EKF, especially in three dimensional and in terrain mapping cases using synchronous or 


asynchronous bearings. 
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APPENDIX A: QWIP SIMULATION PROGRAMS 


In the present Appendix we present selected programs of the simulation of the 
InGaAs/ AlGaAs quantum well. The core of the simulation is E InGaAs.m. Once this 


program is executed, every other sub-matlab file is called automatically. 


%**Simulation Of Performance Of Infrared Photodetectors******* 
% E_ InGaAs 
%LT Psarakis Eftychios Hellenic Navy 
cle 
clear all 
global EE c 
“defining the constants 
hbar = 1.05557e-34; “%[J*s] 
m0 = 9.10939e-31; %[kg] 
q = 1.60218e-19; %[C] 
F = input('Electric Field in [V/m]: "\; 
% Percentage in region 
x=[0,0.3,0.1,0]; 
% Width of the region 
width=[inf,40,40,inf]; 
ii=4; 
for tt=1:11 

%X(i)=input('Percentage of Al (0-1) in region '); 

%width(i) = input('Give width: '); 
end 
y=0; %dlnitial value of y & z to correct the logic of the later condition 
z=0; 
for tt=2:3 
L(1) = width(1)*1e-10; 
L(2) = width(2)* 1e-10; 
L(tt+1) = (width(tt)+width(tt+1))*1e-10; 
end 
if F>-2*1045 & F<10%5 

%EwellNoBiasAiry 

EwellsmallBiasAiry 

%EwellNoBias 

break 
end 
if F>=1045 

for tt=1:4 

if width(tt)==inf 
Z_dir(tt)=0; 
else 
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Z_dir(tt) = L(tt); 
end 
end 
Z_ dir; 
end 
if F<=-2* 1045 
for tt=1:4 
if width(tt)==inf 
Z_dir(tt)=0; 
else 
Z_dir(1) = L(3); 
Z_dir(2) = L(2); 
end 
end 
x=[x(4),x(3),x(2),x(1)] 
end 
Z_dir 
eg(1) = 1.519 + 1.247.*x(1); 
eg(2) = 1.519 - 1.102.*x(2); 
eg(3) = 1.519 - 1.102.*x(3); 
eg(4) = 1.519 + 1.247.*x(4); 
eg=[eg(1) eg(2) eg(3) eg(4)]; 
for tt=1:1 
if F>0 
V(tt) = 0.62*(eg(tt) - eg(2)); 
else 
V(tt) = 0.62*(eg(tt) - eg(3)); 
end 
m_eff(tt) = 1/((x(tt)/(0.028*m0)) + ((1-x(tt))/(0.067*m0))); 
C(tt) = ((2*m_eff(tt)*q)/((F*hbar)*2))(1/3); 
end 
for tt=1:11-1 
sigma(tt) = ((m_eff(tt)*C(tt+1))/(m_eff(tt+1)*C(tt))); 
end 
if F <0; 
DD = 0; 
V2_=V(i); 
else 
DD = -F*L (2); 
V2_ = V(1)- F*L@); 
end 
WEEEEEEEEEEEEEEEEKKE&&continous statesKLKKEEKLKEKEEEKEEE 
ac =1; 
n =0; 
li=1; 
%by varing dE we can ajust the number of states found 
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E_c =min(V(1),V(4))-.03; 
Vb=input('give Vb in eV :'); 
disp(‘The continuus states energy are:') 
% dE_c = 1.2490e-004; 
dE c=(Vb-E_c)/1001; 
while E_c< Vb 
E c=E c+dE ¢c; 
for jj=1:4 
eta(jj) = (E_c - VQj)); 
end 
for jj=1:4 
Alpha(jj) = -CGji)*(F*Z_dir(jj) + eta(jj)); 
Alphaout(11,jj)=[Alpha(:,jj)]; 
end 
for jj=1:3 
alpha(1)=0; 
alpha(Qjj+1)= -CGj+1)*(F*Z_dirGj) + etaGj+1)); 
alphaout(II,jj+1)=[alpha(:,jj+1)]; 
end 
for jj=1:3 
ml11_c=(pi)*(airy(0,alpha(Gjj+1))*airy(3,Alpha(jj))- 
sigma(jj)*(airy(2,Alpha(jj))*airy(1,alpha(jj+1)))); 
ml12_c= (pi)*(airy(2,alpha(jj+1))*airy(3,Alpha(jj)) - 
sigma(jj)*(airy(2,Alpha(jj))*airy(3,alpha(jj+1)))); 
m21_c= (pi)*(sigma(j)*(airy(0,Alpha(jj))*airy(1,alpha(j+1))) : 
airy(1,Alpha(jj))*airy(0,alphaQj+1))); 
m22_c= (pi)*(sigma(j)*(airy(0,Alpha(jj))*airy(3,alpha(j+1))) : 
airy(1,Alpha(jj))*airy(2,alphaQj+1))); 
“checking the output data 
ml lout_c(1L,jj)=[m11_c]; 
m1 2out_c(IL,jj)=[m12_c]; 
m2 lout_c(IL,jj)=[m21_c]; 
m22out_c(Il,jj)=[m22_c]; 
M_c(:jj)=[m11_c;m12_c3;m21_c;m22_c]; 
end 
“%ofinding Mn 
M12_c=[M_c(1,1),M_c(2,1);M_c(3,1),M_c(4,1)]*[M_c(1,2),M_c(2,2);M_c(3,2), 
M_c(4,2)]; 
M23 _c=[M_c(1,2),M_c(2,2);M_c(3,2),M_c(4,2)]*[M_c(1,3),M_c(2,3);M_c(3,3), 
M_c(4,3)]; 
M123_c=[M_c(1,1),M_c(2,1);M_c(3,1),M_c(4,1)]*[M_c(1,2),M_c(2,2);M__c(3,2) 
»M_c(4,2)]*[M_c(1,3),M_c(2,3);M_c(3,3),M_c(4,3)]; 
disp(E_c) 
Eout_cont(ac) = E_c; 
Epsi_cont; 
hold on 
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r = 0:.1:width(2)+width(3); 
plot(r,Eout_cont(ac),'r') 
text(0,Eout_cont(ac)+.01,['E_';num2str(ac+2),' 
‘"sprintf('%4.3f,Eout_cont(ac)),' eV']) 
plot((X*1e10),(real(Psil_cont(ac,:)/5e5)+E_c),'k') 
plot((Y*1e10),(real(Psi2_cont(ac,:)/5e5)+E_c),'k') 
plot((Z*1e10),(real(Psi3_cont(ac,:)/5e5)+E_c),'k') 
plot((T*1e10),(real(Psi4_cont(ac,:)/5e5)+E_c),'k') 
ac=act1; 
ll=I1+1; 
end 
text(-40,1.15*V(1),['V(1) = ',sprintf('%4.3f,V(1)),' eV"]) 
xlabel("X (Angstroms)') 
ylabel('Potential Height (eV)') 
title((Quantum Well’) 
WEEEEEEEEEEEEEEKK&&Endofcontinuous statesKLKKKEKEEEE 
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hold on 

xl = -100:.1:0; 

hl = 0:.001:V(1); 

plot(x1,(V(1)-F*(x1* le-10)),'LineWidth',3) 


x2 = 0:.1:width(2); 
x3 = width(2):.1:(width(2)+width(3)); 


ifF>0 
h2 = (-F*L(2)):.001:(V(3)-(F*L(2))); 
h3 = (V(3)-(F*L@3))):.001:(V(1)-(F*L(3))); 
plot(x3,((V(3)-F*(x3*1e-10))) ,"LineWidth',3) 
end 


ifF <0 

h2 = (-F*L(2)):.001:(V(2)-F*L(2)); 

h3 = (V(2)-F*L(3)):.001:(V(1)-F*L(3)); 

plot(x3,((V(2)-F*(x3*1le-10))) ,"LineWidth',3) 
end 
x4 = (width(2)+width(3)):.1:(width(2)+width(3)+100); 
plot(x4,((V(1)-F*(x4* 1le-10))),'Line Width',3) 
plot(x2,(-F*(x2*1e-10)),'LineWidth',3) 
plot(0,h1,'Line Width', 1) 
plot(width(2),h2,'LineWidth',1) 
plot((width(2)+width(3)),h3,'Line Width’, 1) 
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% 1/7 ee 2 2 2 6 6 26 2 2 6 6 6 2 2 2 6 6 262 2 6 6 Ae 2 2 ae fe 6 2 2 2 fe he Ae 2 2 ae fe 6 2 2 2g fe ie 2 2 2s ie fe ie 2 2k 2k 
% 1/7 He ie 8 2 2 2 ae he 26 2 2 ae fe 6 2 2 2 6 6 26 2 2 6 26 Ae 2 2 2g 6 6 2 2 2g 6 he Ae 2 2 2 fe 6 2 2 2g fe ie 2 2 2s ae fe ie 2 2 2k 


clear an Il dE E eta alpha Alpha E; 

a=l; 

n =0; 

li=1; 

dE = 5*10*-4; 

E =DD; 

%Evaluate the energy from the bottom to the top of the well; 
while E < V2_-dE 

if E>V2_-10%-3 











dE = 10*-6; 

end 

E=E+dE; 

for tt=1:4 

eta(tt) = (E - V(tt)); 

end 

for jj=1:4 

Alpha(jj) = -C@ji)*(F*Z_dir(jj) + eta(jj)); 

end 

for jj=1:3 

alpha(jj+1)= -CGj+1)*(F*Z_dirGj) + etaGj+1)); 

end 

for jj=1:3 

ml1= (pi)*(airy(0,alpha(jj+1))*airy(3,Alpha(jj)) 
sigma(jj)*(airy(2,Alpha(jj))*airy(alpha(jj+1)))); 

m12= (pi)*(airy(2,alpha(jj+1))*airy(3,Alpha(jj)) 
sigma(jj)*(airy(2,Alpha(jj))*airy(3,alpha(jj+1)))); 

m21= (pi)*(sigma(ij)*(airy(0,Alpha(jj))*airy(1,alphaGj+1))) 
airy(1,AlphaGj))*airy(0,alpha(jj+1))); 

m22= (pi)*(sigma(ij)*(airy(0,Alpha(jj))*airy(3,alphaGj+1))) 
airy(1,AlphaGj))*airy(2,alpha(jj+1))); 





“checking the output data 

m1 lout(IL,jj)=[m1 1]; 

m1 2out(IL,jj)=[m12]; 

m2 1out(I1,jj)=[m2 1]; 

m22out(Il,jj)=[m22]; 

M(:,jj)=[m11;m12;m21;m22]; 

end 

“finding Mn 

M12=[M(1,1),M(2,1);M(3,1),M(4,1)]*[M(1,2),M(2,2);M(3,2),M(4,2)]; 

M23=[M(1,2),M(2,2);M(3,2),M(4,2)]*[M(1,3),M(2,3);M(3,3),M(4,3)]; 

M123=[M(1,1),M(2,1);M(3,1),M(4,1)]*[M(1,2),M(2,2);M(3,2),M(4,2)]*[M(1,3), 
M(2,3);M(3,3),M(4,3)]; 
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M22 = real(M123(2,2)); %in case we have three matrices 
“checking the output data 
M22out(Il)=[M22]; 
EM22(11)=[E]; 
MWWW%W%%Y%EN of check%%%%%%%%%%%%%%% 
z= Yy; 
y = M22; 
ll=11+1; 
if z*y <0 “Determine the approximate eigen energy first 
n=n+l; 
ifn<2 “Increase the acuaracy of eigen energy 
E = E-dE; 
dE = le-6; 
y=Z, 
end 
end 
ifn>1 
disp(‘The ground state energy is:') 
disp(E) 
Eout(a) = E; Stores energy levels into an array 
Epsi_bound; 
dE = 5e-4; 
n= 0; 
hold on 
r = 0:.1:width(2)+(a-1)*width(3); 
ifa<3 
plot(r,Eout(a),'r') 
text(0,Eout(a)+.01,['E_'snum2str(a),' = ',sprintf('%4.3f,Eout(a)),' eV']) 
plot((X*1e10),(real(Psil_ bound(a,:)/S5e5)+Eout(a)),'k') 
plot((Y*1e10),(real(Psi2_bound(a,:)/S5e5)+Eout(a)),'k') 
plot((Z*1e10),(real(Psi3_bound(a,:)/5e5)+Eout(a)),'k') 
plot((T*1e10),(real(Psi4_bound(a,:)/5e5)+Eout(a)),'k') 
end 
a=atl; 
end 
end 
ZETA; 


1/7 2 He ae he He 2 2 2 6 6 2 2 2 6 he 26 2 2 ee 6 2 2 2 6 6 26 2 2 ae 6 Ae 2 2 2 6 6 2 2 2 6 246 Ae 2 2g ee ie 2c 2s 2k fe ie Ie 2k 








1/7 2s He ae He 2 2 2 EEK GQ Oi] ator Streriath *F#+Fttett ter eeee eee eee 
1/7 2 2g ie he he 2 2 2g fe 8 2 2 2 fe 6 26 2 2 ae fe 6 2 2 2 fe 6 26 2 2 6 246 26 2 2 2g 6 6 2 2 2 6 fe Ae 2 2s ae fe ie 2s 2s 2k 6 ie 2 2k 
for ai=1:ac-1; 
DE(ai)=(Eout_cont(ai)-Eout); 
f oscil(ai)=2*(m_eff(1))/(hbar’2).*DE(ai)*q.*Zeta_final(ai)'; 
end 
figure 
plot(DE,real(f_oscil)) 
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xlabel(‘\DeltaE') 
ylabel(‘oscillator strength f') 
title(‘oscillator strength Vs \DeltaE') 

for ai=1:ac-1; 
Prod _f oscil _ DOS(ai)=f_oscil(ai).*DOS(ai)*q; 
OSCSTR(ai)=f_oscil(ai).*DOS(ai)*q*dE_c; 
OSCSTR_INTEG=sum(OSCSTR); 

end 
figure 
plot(DE,real(Prod_f oscil DOS)) 
xlabel(‘\DeltaE') 
ylabel(‘(oscillator strength)*(density of states)') 
title(‘oscillator strength multiplied by density of states Vs \DeltaE') 
disp('The (oscillator strength)*(density of states) integral is') 
round(OSCSTR_INTEG) 


% Vent rer ent eee aDSUI UO = te ee er Te Te 

% 1/78 Fe he 8 2 2 2 ae 6 26 2 2 6 fe 6 2 2 2 6 6 2 2 2 6 246 fe ie 26 2 2 6 26 Ae 2 2 ae fe 6 2k 2 2g 6 6 26 2 2 he fe ie 2 2s 2g fe ie 2 2 

q = 1.60218e-19; %[C] 

h=6.626* 10%-34; 

m_eff(1)=1.0e-031*0.6103; 

N_d=1e24; 

e0=8.85e-12; 

n_r=3.5; 

c=3* 1048; 

Cb_c=pi*(N_d*q*hbar’2)./(2*m_eff(1).42*e0*n_r*c); 

total_ length=380*10*-10; 

for kk=1:ac-1 

beta4=((hbar’2)/(2*m_ eff(4)*F))*(1/3); 

“region IV 

quantum _eff=total_length*Cb_c.*Prod_f oscil _DOS; 

J4R(kk)=hbar/(4*m_eff(4)*pi*abs(beta4))*(b3_ cont(kk).*A4(kk)+b3_cont(kk).*( 
i*B4(kk))).*conj((b3_ cont(kk).*A4(kk)+b3_cont(kk).*(i*B4(kk))))*q; 

% I=J_total* Area; 

end 

J_total=quantum_eff.*J4R'.*q; 

figure 

plot((h*c)./(q*DE),real(J_ total)/max(real(J_total)),'b') 

hold on 

plot((h*c)./(q*DE),real(Cb_c*Prod_f oscil _ DOS)./max(real(Cb_c*Prod_f_ oscil_ 
DOS)),'r') 

xlabel("'Wavelength-\mum') 

ylabel('Normalized photocurrent, Normalized absorption’) 

legend('photocurrent','absorption') 
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%**Simulation Of Performance Of Infrared Photodetectors******* 
% EwellNoBias 
%LT Psarakis Eftychios Hellenic Navy 
hbar = 1.05557e-34; “%[J*s] 
m0 = 9.10939e-31; %[kg] 
q = 1.60218e-19; %[C] 
F =F; 
% Percentage in region 
x=[0,0.3,0.1,0]; 
% Width of the region 
width=[inf,40,40,inf]; 
ii=4; 
yao 
z= 0; 
for tt=2:11-1 
L(1) = width(1)* 1le-10; 
L(2) = width(2)* 1e-10; 
L(tt+1) = (width(tt)+width(tt+1))* 1e-10; 
end 
eg(1) = 1.519 + 1.247.*x(1); 
eg(2) = 1.519 - 1.102.*x(2); 
eg(3) = 1.519 - 1.102.*x(3); 
eg(4) = 1.519 + 1.247.*x(4); 
eg=[eg(1) eg(2) eg(3) eg(4)]; 
V(1)=0.62*(eg(1)-eg(2)); 
V(2)=0; 
V(3)=0.62*(eg(3)-eg(2)); 
V(4)=0.62*(eg(4)-eg(2)); 
V=[V(1),V(2),V(3),V(4)]; 
for jj=1:4 
m_eff(jj) = 1/((xGj)/(0.028*m0)) + ((1-x(jj))/(0.067*m0))); 

end 
for t=1:11 

if width(t)==inf 

Z_dir(t)=0; 
else 
Z_dir(t) = L(t); 
end 
end 

ac =1; 
n =0; 
li=1; 
%by varing dE we can ajust the number of states found 
E =max(V(1),V(4)); 
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ji=0; 
Vb=input('give Vb in eV :’) 
dE c=10%-4; 
disp(‘The continuus states energy are:') 
while E<Vb 
E=E+dE c; 
for nn=1:4 
k(nn)=sqrt((2*m_ eff(nn)*(E-V(nn))*q)/hbar’2); 
end 
for nn=1:3 
XI(nn)=(k(nn)*m_eff(nn+1))/(k(nn+1)*m_eff(nn)); 
end 
for jj=1:3 
m1 1_c= 1/2*((1+XI(j))*exp(i*(k(jj)-k(ij+1))*Z_dirGj))); 
m12_c= 1/2*((1-XI(j))*exp(-i*(kGj)+k(j+1))*Z_dirGi)); 
m21_c= 1/2*((1-XI(j))*exp(i* (kGj)+k(j+1))*Z_dirGij)))s 
m22_c= 1/2*((1+XIGj))*exp(-i*(k(j)-kGj+1))*Z_dirGi))); 
M_ c(:,jj=[ml1_c3m12_c;m21_c;m22_ cl]; 
m1 lout_c(Il,jj)=[m11_c]; 
m12out_c(Il,jj)=[m12_c]; 
m21lout_c(Il,jj)=[m21_c]; 
m22out_c(Il,jj)=[m22_c]; 
end 
“finding Mn 
M21 c=[M_ c(1,2),M_c(2,2);M_c(3,2),M_c(4,2)]*[M_c(1,1),M_c(2,1);M_c(3,1), 
M_c(4,1)]; 
M32_c=[M_c(1,3),M_c(2,3);M_c(3,3),M_c(4,3)]*[M_c(1,2),M_c(2,2);M_c(3,2), 
M_c(4,2)]; 
M321_c=[M_c(1,3),M_c(2,3);M_c(3,3),M_c(4,3)]*[M_c(1,2),M_c(2,2);M_c(3,2) 
,M_c(4,2)]*[M_c(1,1),M_c(2,1);M_c(3,1),M__c(4,1)]; 











M22 _ c =real(M321_ c(2,2)); 
“checking the output data 
disp(E) 
Eout_cont(ac) = E; 
integration_cont; 
hold on 
r = 0:.1:width(2)+width(3); 
plot(r,Eout_cont(ac),'r') 
text(0,Eou_contt(ac)+.01,['E_';num2str(ac),' = 
‘" sprintf('%4.3f,Eout_cont(ac)),' eV']) 
plot(reg1*10%10,(real(psil_cont(ac)/S5e5)+Eout_cont(ac)),'k') 
plot(reg2*10%10,(real(psi2_cont(ac)/S5e5)+Eout_cont(ac)),'k') 
plot(reg3*10%10,(real(psi3_cont(ac)/S5e5)+Eout_cont(ac)),'k') 
% plot(reg4*10%10,(real(psi4_cont(ac)/S5e5)+Eout_cont(ac)),'k') 
ac=act1; 
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end 

xlabel('X (Angstroms)') 
ylabel('Potential Height (eV)') 
title((Quantum Well’) 


1/7 2s 2s ghee 2 2 2 fe ie 26 2 2 ae he 6 2 2 28 HE HE YT en] |S 6 6 2 2 2 he he 2 2 2 he he 6 2 2 2g fe ie 2 2 2g he fe ie 2 2s 2s fe ie 
1/7 2 He ae he He 2 2 2 6 6 2 2 2 6 6 26 2 2 ae 6 6 2 2 2g 6 6 26 2 2 6 6 Ae 2 2 2 6 6 2 2 2 6 fe Ae 2 2g ae fe ie 2 2 2k 6 ie Ie 2k 


hold on 
x1 = -100:.1:0; 
hl = 0:.001:V(1); 
plot(x1,(V(1)-F*(x1*1e-10))) %,'LineWidth',5) 
x2 = 0:.1:width(2); 
x3 = width(2):.1:(width(2)+width(3)); 
h2 = (-F*L(2)):.001:(V(3)-(F*L(2))); 
h3 = (V(3)-(F*L@G))):.001:(V(1)-(F*L(3))); 
plot(x3,((V(3)-F*(x3*1le-10)))) %,'LineWidth',5) 
x4 = (width(2)+width(3)):.1:(width(2)+width(3)+100); 
plot(x4,((V(1)-F*(x4* 1e-10)))) %,'LineWidth',5) 
x5 = 0:.1:(width(2)+width(3)); 
plot(x2,(-F*(x2*1e-10))) %,'LineWidth',5) 
plot(0,h1) %,'Line Width',8) 
plot(width(2),h2) %,'LineWidth',3) 
plot((width(2)+width(3)),h3) %,'LineWidth',5) 
1/7 2s He ag he he 2 2 2 6 6 2k 2 2 6 46 26 2 2 2 6 6 Fe 2 2 6 he 6 2 2 ee 6 2 2 2 6 46 Ae 2 2 ae fe 6 2s 2 2 fe he 28 2 2s ae fe ie 2 2 2s 6 ie ie 2k 
1/7 2s He ae he he 2 2 2g fe 6 2k 2 2 6 he 26 2 2 ee 6 2 2 2 6 6 Ae 2 2 6 6 6 2 2 2 6 he 26 2 2 ae fe 6 2 2 2g fe 26 26 2 2s 2g fe ie 2s 2s 2k fe ie Ie 2k 
1/57 2 He ae ie He 2 2s 2g he He HE HE AE AED CI States * 7 A A A 2 He He ae He 2 He 2s he ee 2 He ie fe He 2 2k 2k 
a=l; 
n =0; 
li=1; 
dE = 5*10%-4; 
E =0; 
jj=0; 
while E<max(V) 
if E > max(V)-10%-3 
dE = 10%-6; 
end 
E=E+dE; 
for nn=1:4 
k(nn)=sqrt((2*m_ eff(nn)*(E-V(nn))*q)/hbar’2); 
end 
for nn=1:3 
XI(nn)=(k(nn)*m_ eff(nn+1))/(k(nn+1)*m_ eff(nn)); 
end 
for jj=1:3 
m1 1= 1/2*((1+X1(j))*exp(i*(kGj)-kGi+1))*Z_dirGi)) 
m12= 1/2*((1-XIGj))*exp(-i*(k(ji)+kGi+))*Z_dir(jj))); 
m21= 1/2*((1-X1Gj))*exp(i*(k(i)+kGj+1))*Z_dirGi))); 
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m22= 1/2*((1+XIGj))*exp(-i*(k(i)-kGj+1))*Z_dir(ij))); 
M(:,jj)=[m11;m12;m21;m22]; 
m1 lout(IL,jj)=[m1 1]; 
m1 2out(IL,jj)=[m12]; 
m2 1out(I1,jj)=[m2 1]; 
m22out(I1,jj)=[m22]; 
end 
“finding Mn 
M21=[M(1,2),M(2,2);M(3,2),M(4,2)]*[M(1,1),M(2,1);M(3,1),M(4,1)]; 
M32=[M(1,3),M(2,3);M(3,3),M(4,3)]*[M(1,2),M(2,2);M(3,2),M(4,2)]; 
M321=[M(1,3),M(2,3);M(3,3),M(4,3)]*[M(1,2),M(2,2);M(3,2),M(4,2)]*[M(1,1), 
M(2,1);M(3,1),M(4,1)]; 
M22 =real(M321(2,2)); 
“checking the output data 
M22out(1l)=[M22]; 
EM22(11l)=[E]; 
MWWV%%%%EN of check%%%%%%%%%%%%%%% 
ZY; 
y=M2? ; 
lI=11+1; 
if z.*y <0 
n=nt+l; 
ifn<2 “Increase the acuaracy of eigen energy 
E = E-dE; 
dE = le-6; 
y=Z; 
end 
end 
ifn>1 
disp('The ground state energy is:') 
disp(E) 
Eout(a) = E; 
Xx=E; 
integration; 
dE = 5e-4; 
n= 0; 
hold on 
r = 0:.1:width(2)+(a-1)*width(3); 
Eout(a) = E; 
hold on 
if a<3 
plot(r,Eout(a),'r') 
text(0,Eout(a)+.01,['E 'snum2str(a),' = ',sprintf('%4.3f,Eout(a)),' eV']) 
plot(reg1*1e10,(real(Psil_ bound/5e5)+E),'k') 
plot(reg2*1e10,(real(Psi2_bound/5e5)+E),'k') 
plot(reg3*1e10,(real(Psi3_ bound/5e5)+E),'k') 
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plot(reg4*1e10,(real(Psi4_bound/5e5)+E),'k') 
end 
a=atl; 
end 
end 
ZETAIL: 


Oy 2 Ee 2s 6 2 2 6 2 2 6 2 2 fe 2A 2 6 2 2 fe Ae 2 6 2 2 6 2 2 6 2 2 6 2 2 6 2 2 fe 2A 2 fe 2A 2 fe 2A 2 fe 2 2 fe 2 2 fe 2 2 ie 
Oy 28 2 2 He He 2 2 2 2 EE EK Cy] lator Strength? tte Ft sete settee tes 
1/7 2s 2s ae hee 2 2 2g fe 8 2 2 2 6 6 26 2 2 ae fe he 2 2 2 fe 6 2 2 2 6 46 26 2 2 ae fe 6 2 2 2 6 ie 2 2 2s fe fe ie 2 2s ae fe ie 2 2k 
for ai=1:ac-1; 
DE(ai)=(Eout_cont(ai)-Eout(1)); 
f oscil(ai)=2*(m_eff(1))/(hbar’2).*DE(ai)*q.*Zeta_final(ai)'; 
end 
figure(2) 
plot(DE,real(f_oscil)) 


xlabel(‘\DeltaE') 
ylabel(‘oscillator strength f') 
title(‘oscillator strength Vs \DeltaE') 

for ai=1:ac-1; 
Prod f oscil _ DOS(ai)=f_oscil(ai).*DOS(ai)*q; 
OSCSTR(ai)=f_oscil(ai).*DOS(ai)*q*dE_c; 
OSCSTR_INTEG=sum(OSCSTR); 

end 
figure(3) 
plot(DE,real(Prod_f oscil DOS)) 
xlabel(‘\DeltaE') 
ylabel(‘(oscillator strength)*(density of states)') 
title(‘oscillator strength multiplied by density of states Vs \DeltaE') 
disp('The (oscillator strength)*(density of states) integral is') 
round(OSCSTR_INTEG) 


%**Simulation Of Performance Of Infrared Photodetectors******* 
% integration 
%LT Psarakis Eftychios Hellenic Navy 
dx = 10*-11; 
reg1= -300e-10:dx:0; 
reg2 = 0:dx:L(1,2); 
reg3 = L(1,2):dx:L(1,3); 
reg4 = L(1,3):dx:(L(1,3)+300e-10); 


Al1=0; 
Bl=1; 
A2=M(2,1); 
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B2=M(4,1); 

A3=M21(1,2); 

B3=M21(2,2); 

A4=M321(1,2); 

B4=0; 
psil_bound=A1*exp(i*k(1)*reg1)+B1*exp(-i*k(1)*reg1); 
psi2_bound=A2*exp(i*k(2)*reg2)+B2*exp(-i*k(2)*reg2); 
psi3_ bound=A3*exp(i*k(3)*reg3)+B3*exp(-i*k(3)*reg3); 
psi4_ bound=A4*exp(i*k(4)*reg4)+B4*exp(-i*k(4)*reg4); 
integrall=sum((psil_bound.*conj(psil_bound))*dx); 
integral2=sum((psi2_bound.*conj(psi2_bound))*dx); 
integral3=sum((psi3_bound.*conj(psi3_ bound))*dx); 
integral4=sum((psi4_bound.*conj(psi4_bound))*dx); 
total_integral=integral1+integral2+integral3+integral4; 
b1=(sqrt(total_integral)); 

Psil_bound=1/b1*psil_ bound; 

Psi2_bound=1/b1*psi2_ bound; 

Psi3_bound=1/b1*psi3_ bound; 
Psi4_bound=1/b1*psi4_bound; 

Integrall=sum((Psil_ bound.*conj(Psil_bound))*dx); 
Integral2=sum((Psi2_bound.*conj(Psi2_bound))*dx); 
Integral3=sum((Psi3_bound.*conj(Psi3_bound))*dx); 
Integral4=sum((Psi4_bound.*conj(Psi4_bound))*dx); 
Total_integral=Integral1+Integral2+Integral3+Integral4; 





%**Simulation Of Performance Of Infrared Photodetectors******* 
% integration_cont 
%LT Psarakis Eftychios Hellenic Navy 
dx = 10-11; 

reg1= -300e-10:dx:0; 

reg2 = 0:dx:L(1,2); 

reg3 = L(1,2):dx:L(1,3); 

reg4 = L(1,3):dx:(L(1,3)+300e-10); 
Al_c(ac)=1/sqrt(2*pi); 
B1_ c(ac)=-Al_c(ac)*M321_c(2,1)/M321_ c(2,2); 
A2_c(ac)=M_c(1,1)*Al_c(ac)+M_c(2,1)*B1_ c(ac); 
B2_c(ac)=M_c(3,1)*Al_c(ac)+M_c(4,1)*B1_c(ac); 
A3_c(ac)=M21_ c(1,1)*A1_ c(ac)+M21_ c(1,2)*B1_ c(ac); 
B3_c(ac)=M21_ c(2,1)*Al1_ c(ac)+M21_ c(2,2)*B1_ c(ac); 
A4 c(ac)=M321 c(1,1)*Al_c(ac)+M321 c(1,2)*B1_ c(ac); 
B4 c(ac)=0; 
psil_cont(ac,:)=A1_c(ac)*exp(i*k(1)*reg1)+B1_c(ac)*exp(-i*k(1)*reg1); 
psi2_cont(ac,:)=A2_c(ac)*exp(i*k(2)*reg2)+B2_c(ac)*exp(-i*k(2)*reg2); 
psi3_cont(ac,:)=A3_c(ac)*exp(i*k(3)*reg3)+B3_c(ac)*exp(-i*k(3)*reg3); 
psi4_cont(ac,:)=A4_c(ac)*exp(i*k(4)*reg4)+B4_c(ac)*exp(-i*k(4)*reg4); 
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%**Simulation Of Performance Of Infrared Photodetectors******* 
% ZETAIL 
%LT Psarakis Eftychios Hellenic Navy 
dx = 10%-11; 

reg1= -300e-10:dx:0; 

reg2 = 0:dx:L(1,2); 

reg3 = L(1,2):dx:L(1,3); 

reg4 = L(1,3):dx:(L(1,3)+300e-10); 

dim=max(reg4)-min(reg 1); 

for kk=1:ac-1 
Zital(kk,:)=conj(Psil_bound(1,:)).*reg1.*psil_cont(kk,:).*dx; 
Zita2(kk,:)=conj(Psi2_bound(1,:)).*reg2.*psi2_cont(kk,:).*dx; 
Zita3(kk,:)=conj(Psi3_bound(1,:)).*reg3.*psi3_cont(kk,:).*dx; 
Zita4(kk,:)=conj(Psi4_bound(1,:)).*reg4.*psi4_cont(kk,:).*dx; 
OSC_integral1(kk,:) = sum(Zital(kk,:)); 
OSC_integral2(kk,:) = sum(Zita2(kk,:)); 
OSC_integral3(kk,:) = sum(Zita3(kk,:)); 
OSC_integral4(kk,:) = sum(Zita4(kk,:)); 
OSC_finali(kk,:)}=(OSC_ integral 1(kk,:)+ OSC_integral2(kk,:)+ 
OSC_integral3(kk,:)+OSC_integral4(kk,:)); 
Zeta_final(kk,:)=(OSC_finali(kk,:)).*conj(OSC_finali(kk,:)); 
%DOS(kk)=dim./(hbar*pi)*sqrt(m_eff(1)./(2*q*(Eout_cont(kk)- 
Eout_cont(1)+0.001))); 
DOS(kk)=1/(hbar*pi)*sqrt(m_eff(1)./(2*q*(Eout_cont(kk)-V(1)+0.0001))); 
end 


%**Simulation Of Performance Of Infrared Photodetectors******* 
%EwellsmallBiasAiry 
%LT Psarakis Eftychios Hellenic Navy 
if F>0 
for tt=1:4 
if width(tt)==inf 
Z_dir(tt)=0; 
else 
Z_dir(tt) = L(tt); 
end 
end 
Z_ dir 
else 
for tt=1:4 
if width(tt)==inf 
Z_dir(tt)=0; 
else 
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Z_dir(1) = L(3); 
Z_dir(2) = L(2); 
end 
end 
x=[x(4),x(3),x(2),x(1)]; 


end 


eg(1) = 1.519 + 1.247.*x(1); 
eg(2) = 1.519 - 1.102.*x(2); 
eg(3) = 1.519 - 1.102.*x(3); 
eg(4) = 1.519 + 1.247.*x(4); 
eg=[eg(1) eg(2) eg(3) eg(4)]; 
for tt=1:11 
if F>0 
V(tt) = 0.62*(eg(tt) - eg(2)); 
else 
V(tt) = 0.62*(eg(tt) - eg(3)); 
end 
m_eff(tt) = 1/((x(tt)/(0.028*m0)) + ((1-x(tt))/(0.067*m0))); 
C(tt) = (((2*m_eff(tt)*q)/((F*hbar)*2))*(1/3)); 
end 
for tt=1:11-1 
sigma(tt) = ((m_eff(tt)*C(tt+1))/(m_eff(tt+1)*C(tt))); 
end 


if F <0; 

DD = 0; 
vV2_=V(i); 

else 

DD = -F*L(2); 

V2_ = V(1)- F*L@); 
end 

ac =1; 

n = c 

ll=1; 


%by varing dE we can ajust the number of states found 
E_c =min(V(1),V(4))-.03; 
Vb=input('give Vb in eV :'); 

disp('The continuus states energy are:') 
dE c=(Vb-E _c)/1001; 

while E_c< Vb 

E c=E c+dE ¢c; 

for tt=1:4 

eta(tt) = (E_c - V(tt)); 

end 

for jj=1:4 

Alpha(jj) = -C(jj)*(F*Z_dir(jj) + eta(jj)); 
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Alphaout(I1,jj)=[Alpha(:.jj)]; 

end 

for jj=1:3 

alpha(jj+1)= -CGj+1)*(F*Z_dir(jj) + etaGj+1)); 
alphaout(11,jj)=[alpha(: jj) ]; 











end 

for jj=1:3 

mll c= (pi)*(airy(0,alpha(jj+1))*airy(3,Alpha(jj)) - 
sigma(jj)*(airy(2,Alpha(jj))*airy(1,alpha(jj+1)))); 

ml2_ c= (pi)*(airy(2,alpha(Qj+1))*airy(3,Alpha(jj)) - 
sigma(jj)*(airy(2,Alpha(jj))*airy(3,alpha(jj+1)))); 

m21_ c (pi)*(sigma(jj)*(airy(0,Alpha(jj))*airy(1,alpha(jj+1))) - 
airy(1,Alpha(jj))*airy(0,alphaGj+1))); 

m22 c (pi)*(sigma(jj)*(airy(0,Alpha(jj))*airy(3,alpha(jj+1))) - 
airy(1,Alpha(jj))*airy(2,alphaQj+1))); 





“checking the output data 
m1 lout_c(Il,jj)=[m11_c]; 
m1 2out_c(Il,jj)=[m12_c]; 
m21lout_c(Il,jj)=[m21_c]; 
m22out_c(Il,jj)=[m22_c]; 
M_c(:,jj=[ml1_e3;m12_c;m21_c3m22_ cl]; 
end 
“finding Mn 
M12 c=[M_ c(1,1),M_c(2,1);M_c(3,1),M_c(4,1)]*[M_c(1,2),M_c(2,2);M_c(3,2), 
M_c(4,2)]; 
M23 _c=[M_c(1,2),M_c(2,2);M_c(3,2),M_c(4,2)]*[M_c(1,3),M_c(2,3);M_c(3,3), 
M_c(4,3)]; 
M123 c=[M_ c(1,1),M_c(2,1);M_c(3,1),M_c(4,1)]*[M_c(1,2),M_c(2,2);M_c(3,2) 
,M_c(4,2)]*[M_c(1,3),M_c(2,3);M_c(3,3),M_c(4,3)]; 
disp(E_c) 
Eout_cont(ac) = E_c; 
Epsi_cont; 
hold on 
r = 0:.1:width(2)+width(3); 
plot(r,Eout_cont(ac),'r') 
text(0,Eout_cont(ac)+.01,['E_'jnum2str(ac+2),' = 
‘"sprintf('%4.3f,Eout_cont(ac)),' eV']) 
plot((X*1e10),(real(psil_ cont(ac,:)/5e5)+Eout_cont(ac)),'k') 
plot((Y*1e10),(real(psi2_cont(ac,:)/5e5)+Eout_cont(ac)),'k') 
plot((Z*1e10),(real(psi3_cont(ac,:)/5e5)+Eout_cont(ac)),'k') 
plot((T*1e10),(real(psi4_cont(ac,:)/5e5)+Eout_cont(ac)),'k') 
ac=act1; 
I=11+1; 
end 
text(-40,1.15*V(1),['V(1) = "sprintf('%4.3f,V(1)),' eV']) 
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xlabel('X (Angstroms)') 
ylabel('Potential Height (eV)') 
title((Quantum Well’) 


Of 0/578 2 ae He 2 2s 2s ae he 26 2 He ae fe He 2 2 2g 6 6 AS IF yea] |) 26 6 A 2 2g he he 6 2 2 2g he he 6 2 2 he fe 26 2 2 2s ae 26 ie 2 2s 28 
Of 0/578 A ae He 2 2 2g 6 6 Ae 2 2 2g he 6 2 2 2 6 6 26 2 2 ae he 6 2 2 2g 6 6 2A 2 2 6 he Ae 2 2 ae fe 6 2 2 2g fe ie 2 2 2s ie he ie 2 2 2k 


hold on 

x1 = -100:.1:0; 

hl = 0:.001:V(1); 

plot(x1,(V(1)-F*(x1*1e-10)),'LineWidth',3) 

x2 = 0:.1:width(2); 

x3 = width(2):.1:(width(2)+width(3)); 

ifF >0 

h2 = (-F*L(2)):.001:(V(3)-(F*L(2))); 

h3 = (V(3)-(F*L@G))):.001:(V(1)-(F*L(3))); 

plot(x3,((V(3)-F*(x3*1e-10))) ,"LineWidth',3) 

end 

ifF <0 

h2 = (-F*L(2)):.001:(V(2)-F*L(2)); 

h3 = (V(2)-F*L(3)):.001:(V(1)-F*L(3)); 

plot(x3,((V(2)-F*(x3*1e-10))) ,"LineWidth',3) 

end 

x4 = (width(2)+width(3)):.1:(width(2)+width(3)+100); 

plot(x4,((V(1)-F*(x4* 1e-10))),"LineWidth',3) 
plot(x2,(-F*(x2*1e-10)),'LineWidth',3) 

plot(0,h1,"Line Width’, 1) 

plot(width(2),h2,'LineWidth',1) 

plot((width(2)+width(3)),h3,'Line Width’, 1) 

Oo) 7 He He 2 2 2 3s he 26 2 2 2g fe 6 2 2 2 6 6 26 2 2 2 fe 6 2g 2 2 6 46 Ae 2 2 ae 6 6 2 2 2 6 246 26 2 2g ee 6 2 2 2 6 ie 26 2 2s 6 26 ie 2 2s 2s Ie ik 
1/7 2s He ae ie he 2 2s ae he He 2 2 2 EEE BOUIN DD) ST AT ES 5 8 28 28 28 26 2 2 2 2s 2s he 2 2 2g ae ee 2 2s 2c fe ie 2 2s 2s Ie ie 
1/7 2 He ae he 26 2 2g 2g 6 6 2 2 2 6 he 26 2 2 2 fe 6 Ze 2 2 6 6 Ae 2 2 ee 6 2 2 2 6 6 Ae 2 2 2 fe 6 22 2 fe he 26 2 2 2c fe ie 2 2 2k fe ie 2s 2k 


% clear M22 M22out n Il dE E m11 m12 m21 m22 alpha alphaout 
% clear m1 lout m12out m21out m220ut M12 M23 M123 Alpha Alphaout 
ab =1; 


n =0; 

ll=1; 

dE = 5*10%-4; 
E =DD; 


%Evaluate the energy from the bottom to the top of the well; 
while E < V2_-dE 


if E> V2_-10%-3 
dE = 10%-6; 

end 

E=E+dE; 

for tt=1:4 

eta(tt) = (E - V(tt)); 
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end 
for jj=1:4 
Alpha(jj) = -CGji)*(F*Z_dir(jj) + eta(jj)); 
Alphaout(11,jj)=[Alpha(:,jj)]; 
end 
for jj=1:3 
alpha(jj+1)= -CGj+1)*(F*Z_dirGj) + etaGj+1)); 
alphaout(I1,jj)=[alpha(:,jj)]; 
end 
for jj=1:3 
m1 1=(pi)*(asymptoticAiry(0,alpha(jj+1))*asymptoticAiry(3,Alpha(jj)) - 
sigma(jj)*(asymptoticAiry(2,Alpha(jj))*asymptoticAiry(1,alpha(jj+1)))); 
m12=(pi)*(asymptoticAiry(2,alpha(jj+1))*asymptoticAiry(3,Alpha(jj)) - 
sigma(jj)*(asymptoticAiry(2,Alpha(jj))*asymptoticAiry(3,alpha(jj+1)))); 
m21= 
(pi)*(sigma(jj)*(asymptoticAiry(0,Alpha(jj))*asymptoticAiry(1,alpha(Gj+1))) - asymp- 
toticAiry(1,Alpha(jj))*asymptoticAiry(0,alpha(jj+1))); 
m22= 
(pi)*(sigma(jj)*(asymptoticAiry(0,Alpha(jj))*asymptoticAiry(3,alpha(Gj+1))) - asymp- 
toticAiry(1,Alpha(jj))*asymptoticAiry(2,alpha(jj+1))); 
“checking the output data 
m1 lout(IL,jj)=[m1 1]; 
m1 2out(IL,jj)=[m12]; 
m2 1out(I1,jj)=[m2 1]; 
m22out(I1,jj)=[m22]; 
M(:,jj)=10%-40*[m1 1;m12;m21;m22]; 
end 
“finding Mn 
M12=[M(1,1),M(2,1);M(3,1),M(4,1)]*[M(1,2),M(2,2);M(3,2),M(4,2)]; %10*60 
M23=[M(1,2),M(2,2);M(3,2),M(4,2)]*[M(1,3),M(2,3);M(3,3),M(4,3)]; %10°60 
M123=[M(1,1),M(2,1);M(3,1),M(4,1)]*[M(1,2),M(2,2);M(3,2),M(4,2)]*[M(1,3), 
M(2,3);M(3,3),.M(4,3)]; %10*90 
M22 = real(M123(2,2)); 
“checking the output data 
M22out(I))=[M22]; 
EM22(1)=[E]; 
MWV%%%%Y%end Of check%%%%%%%%%%%%%%% 
Il=11+1; 
Z~Yy; 
y = M22; 
if z*y <0 “Determine the approximate eigen energy first 
n=n+l; 
ifn<2 “Increase the acuaracy of eigen energy 
E = E-dE; 
dE = le-6; 
y=Z; 
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end 
end 
ifn> 1 
disp('The ground state energy is:') 
disp(E) 
%  disp(M123) 
%  disp(M23) 
Eout(ab)=E; %Stores energy levels into an array 
Epsi_bound1; 
dE = 5e-4; 
n=0; 
hold on 
r = 0:.1:width(2)+(ab-1)*width(3); 
if ab <3 
plot(r,Eout(ab),'r') 
text(0,Eout(ab)+.01,['E ';num2str(ab),' = ',sprintf('%4.3f,Eout(ab)),' eV']) 
plot((X*1e10),(real(Psil_ bound(ab,:)/5e5)+E),'k') 
plot((Y*1e10),(real(Psi2_ bound(ab,:)/5e5)+E),'k') 
plot((Z*1e10),(real(Psi3_bound(ab,:)/5e5)+E),'k') 
plot((T*1e10),(real(Psi4_bound(ab,:)/5e5)+E),'k') 
end 
ab = ab+1; 
end 
end 
ZETA: 


1/7 2 2s ie he Ae 2 2 2g fe 6 2 2 2g 6 6 Ae 2 2 ae fe 6 2 2 2 6 6 26 2 2 6 6 6 2 2 ee 6 2 2 2g ae fe Ae 2 2s 2g fe ie 2 2 2s 6 ie 2 2k 








1/5 7 2 He ae ae A 2 2 2 EEK GQ Oi] ator Strength + teeter eee eter eee eee 
1/7 2 He ag he He 2 2 2g fe 6 2k 2 2 6 he 26 2 2 6 fe 6 2 2 2 6 246 26 2 2 6 fe Ae 2 2 2 fe 6 2 2 2 6 26 6 2 2s 2c fe ie 2 2s 2s ae ie 2s 2k 
for ai=1:ac-1; 
DE(ai)=(Eout_cont(ai)-Eout); 
f oscil(ai)=2*(m_eff(1))/(hbar’2).*DE(ai)*q.*Zeta_final(ai)'; 
end 
figure(2) 
plot(DE,real(f_oscil)) 


xlabel(‘\DeltaE') 
ylabel(‘oscillator strength f') 
title(‘oscillator strength Vs \DeltaE') 
for ai=1:ac-1; 
Prod _f oscil_ DOS(ai)=f_oscil(ai).*DOS(ai)*q; 
OSCSTR(ai)=f_oscil(ai).*DOS(ai)*q*dE_c; 
OSCSTR_INTEG=sum(OSCSTR); 
end 
figure(3) 
plot(DE,real(Prod_f oscil DOS)) 
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xlabel(\DeltaE') 
ylabel(‘(oscillator strength)*(density of states)') 
title(‘oscillator strength multiplied by density of states Vs \DeltaE') 


disp('The (oscillator strength)*(density of states) integral is') 
round(OSCSTR_INTEG) 


%**Simulation Of Performance Of Infrared Photodetectors******* 
%Epsi_bound1 
%LT Psarakis Eftychios Hellenic Navy 
warning off 
deltax = 10-11; 
B4=1; 
A4=0; 
A3(ab) = 10*40*M(2,3)*B4; 
B3(ab) = 10°40*M(4,3)*B4; 
A2(ab) = 10*80*M23(1,2); 
B2(ab) = 10°80*M23(2,2)*B4; 
% if ab==1 
% Al(1)=0; 
% else 
Al1(ab) =10420*M123(1,2); 
% end 
B1=0; 


ifF >0 

X = -300e-10:deltax:0; 

Y = 0:deltax:L(1,2); 

Z = L(,2):deltax:L(1,3); 

T = L(1,3):deltax:(L(1,3)+300e-10); 
else 

X = L(1,3):deltax:(L(1,3)+300e-10); 

Y = L(1,2):deltax:L(1,3); 

Z = 0:deltax:L(1,2); 

T = -300e-10:deltax:0; 
end 
P1_bound(ab,:) = -C(1).*(F.*X + eta(1)); 
P2_bound(ab,:) = -C(2).*(F.*Y + eta(2)); 
P3_bound(ab,:) = -C(3).*(F.*Z + eta(3)); 
P4_bound(ab,:) = -C(4).*(F.*T + eta(4)); 








psil_ bound(ab,:) = (Al(:,ab).*asymptoticAiry(0,P1 bound(ab,:)) 
B1.*asymptoticAiry(2,P1_bound(ab,:))); 
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psi2_bound(ab,:) = (A2(:,ab).*asymptoticAiry(0,P2_bound(ab,:)) a 


B2(:,ab).*asymptoticAiry(2,P2_bound(ab,:))); 


psi3_ bound(ab,:) = (A3(:,ab).*asymptoticAiry(0,P3_bound(ab,:)) = 


B3(:,ab).*asymptoticAiry(2,P3_bound(ab,:))); 


psi4_bound(ab,:) = (A4.*asymptoticAiry(0,P4 bound(ab,:)) a 


B4.*asymptoticAiry(2,P4 bound(ab,:))); 


psil_ bound(ab,:) = 100*psil_ bound(ab,:); 

psi2_bound(ab,:) = 10*-100*psi2_bound(ab,:); 

psi3_ bound(ab,:) = 10“-100*psi3_bound(ab,:); 

psi4_bound(ab,:) = 10“-100*psi4_bound(ab,:); 

integrand! bound(ab,:) = (psil_bound(ab,:).*conj(psil_ bound(ab,:)))*deltax; 
integrand2 bound(ab,:) = (psi2_bound(ab,:).*conj(psi2_bound(ab,:)))*deltax; 
integrand3_bound(ab,:) = (psi3_bound(ab,:).*conj(psi3_bound(ab,:)))*deltax; 
integrand4 bound(ab,:) = (psi4_bound(ab,:).*conj(psi4_bound(ab,:)))*deltax; 
integrall_bound(ab,:) = sum(integrand!1_bound(ab,:)); 

integral2 bound(ab,:) = sum(integrand2_ bound(ab,:)); 

integral3_bound(ab,:) = sum(integrand3_bound(ab,:)); 

integral4 _bound(ab,:) = sum(integrand4 bound(ab,:)); 

finali_bound(ab,:) = integrall bound(ab,:) + integral2 bound(ab,:) +  inte- 


gral3_bound(ab,:) + integral4_bound(ab,:); 


b3_bound(ab,:) = 1./(sqrt(finali_bound(ab,:))); 

Psil_bound(ab,:) = b3_bound(ab,:).*psil_bound(ab,:); 

Psi2_bound(ab,:) = b3_bound(ab,:).*psi2_bound(ab,:); 

Psi3_bound(ab,:) = b3_ bound(ab,:).*psi3_ bound(ab,:); 

Psi4_bound(ab,:) = b3_bound(ab,:).*psi4_bound(ab,:); 
Integrall_bound(ab,:)=sum((Psil_ bound(ab,:).*conj(Psil_bound(ab,:)))*deltax); 
Integral2_bound(ab,:)=sum((Psi2_bound(ab,:).*conj(Psi2_bound(ab,:)))*deltax); 
Integral3_bound(ab,:)=sum((Psi3_ bound(ab,:).*conj(Psi3_ bound(ab,:)))*deltax); 
Integral4_bound(ab,:)=sum((Psi4_bound(ab,:).*conj(Psi4_bound(ab,:)))*deltax); 


Ttal_ integral _bound(ab,:)=Integrall_bound(ab,:)+Integral2_bound(ab,:)+Integral3_boun 
d(ab,:)+Integral4_bound(ab,:) 


%**Simulation Of Performance Of Infrared Photodetectors******* 
%Epsi_contl 
%LT Psarakis Eftychios Hellenic Navy 
deltax = 10-11; 
if F>0 
X = -300e-10:deltax:0; 
Y = 0:deltax:L(1,2); 
Z = L(1,2):deltax:L(1,3); 
T = L(1,3):deltax:(L(1,3)+300e-10); 
AA(ac) = 1/sqrt(2*pi); 
B4(ac) = -(M123_c(2,1)/M123_c(2,2))*A4(ac) ; 
A3(ac) = (M_c(1,3)*A4(ac) + B4(ac)*M_c(2,3)); 
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B3(ac) = (M_c(3,3)*A4(ac) + B4(ac)*M_c(4,3)); 
A2(ac) = (M23_c(1,1)*A4(ac)+B4(ac)*M23_c(1,2)); 
B2(ac) = (M23_c(2,1)*A4(ac)+B4(ac)*M23_c(2,2)); 
Al(ac) = (M123 _c(1,1)*A4(ac) + B4(ac)*M123_c(1,2)); 
Bl(ac) = 0; 

else 
X = L(1,3):deltax:(L(1,3)+300e-10); 
Y = L(1,2):deltax:L(1,3); 
Z = 0:deltax:L(1,2); 
T = -300e-10:deltax:0; 
AA(ac) = 1/sqrt(2*pi); 
B4(ac) = -(M123_c(2,1)/M123_c(2,2))*A4(ac) ; 
A3(ac) = (M_c(1,3)*A4(ac) + B4(ac)*M_c(2,3)); 
B3(ac) = (M_c(3,3)*A4(ac) + B4(ac)*M_c(4,3)); 
A2(ac) = (M23_c(1,1)*A4(ac)+B4(ac)*M23_c(1,2)); 
B2(ac) = (M23_c(2,1)*A4(ac)+B4(ac)*M23_c(2,2)); 
Al(ac) = (M123 _c(1,1)*A4(ac) + B4(ac)*M123_c(1,2)); 
Bl(ac) = 0; 

end 

P1_cont(ac,:) = -C(1).*(F.*X + eta(1)); 

P2_ cont(ac,:) = -C(2).*(F.*Y + eta(2)); 

P3_cont(ac,:) = -C(3).*(F.*Z + eta(3)); 

P4 cont(ac,:) = -C(4).*(F.*T + eta(4)); 

















psil_cont(ac,:) = Al(ac)*asymptoticAiry(0,P1_cont(ac,:)) 
Bl(ac)*asymptoticAiry(2,P1_cont(ac,:)); 

psi2_cont(ac,:) = A2(ac)*asymptoticAiry(0,P2_cont(ac,:)) 
B2(ac)*asymptoticAiry(2,P2_cont(ac,:)); 

psi3_cont(ac,:) = A3(ac)*asymptoticAiry(0,P3_cont(ac,:)) 
B3(ac)*asymptoticAiry(2,P3_cont(ac,:)); 

psi4_cont(ac,:) = AA(ac)*asymptoticAiry(0,P4_ cont(ac,:)) 


B4(ac)*asymptoticAiry(2,P4_ cont(ac,:)); 
integrand! cont(ac,:) = (psil_cont(ac,:).*conj(psil_cont(ac,:)))*deltax; 
integrand2_ cont(ac,:) = (psi2_cont(ac,:).*conj(psi2_cont(ac,:)))*deltax; 
integrand3_cont(ac,:) = (psi3_cont(ac,:).*conj(psi3_cont(ac,:)))*deltax; 
integrand4 cont(ac,:) = (psi4_cont(ac,:).*conj(psi4_cont(ac,:)))*deltax; 
integrall_cont(ac,:) = sum(integrand1_cont(ac,:)); 
integral2_cont(ac,:) = sum(integrand2_cont(ac,:)); 
integral3_cont(ac,:) = sum(integrand3_cont(ac,:)); 
integral4_cont(ac,:) = sum(integrand4 cont(ac,:)); 
finali_cont(ac,:) =  integrall cont(ac,:) + integral2 cont(ac,:) 
gral3_cont(ac,:) + integral4_cont(ac,:); 
b3_cont(ac,:) = 1/(sqrt(finali_cont(ac,:))); 
Psil_cont(ac,:) = b3_cont(ac,:)*psil_cont(ac,:); 
Psi2_cont(ac,:) = b3_ cont(ac,:)*psi2_cont(ac,:); 
Psi3_cont(ac,:) = b3_ cont(ac,:)*psi3_cont(ac,:); 
Psi4_cont(ac,:) = b3_cont(ac,:)*psi4_cont(ac,:); 
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+ 





inte- 


To- 


Integrall_cont(ac,:)}=sum((Psil_cont(ac,:).*conj(Psil_cont(ac,:)))*deltax); 
Integral2_cont(ac,:)=sum((Psi2_cont(ac,:).*conj(Psi2_cont(ac,:)))*deltax); 
Integral3_cont(ac,:)}=sum((Psi3_cont(ac,:).*conj(Psi3_cont(ac,:)))*deltax); 
Integral4_cont(ac,:)=sum((Psi4_cont(ac,:).*conj(Psi4_cont(ac,:)))*deltax); 


tal_integral_cont(ac,:)=Integrall_cont(ac,:)+Integral2_cont(ac,:)+Integral3_cont(ac,:)+Int 
egral4_cont(ac,:); 


%**Simulation Of Performance Of Infrared Photodetectors******* 
%Epsi_bound 
%LT Psarakis Eftychios Hellenic Navy 
deltax = 10-11; 
A4 b=0; 
B4 b= 1; 
A3_b(a) = M(1,3)*A4_b + M(2,3)*B4 _b; 
B3_b(a) = M(3,3)*A4_ b+ M(4,3)*B4 b; 
A2_b(a) = M23(1,2)*B4 b; 
B2_b(a) = M23(2,1)*A4_b + M23(2,2)*B4 b; 
Al_b(a) = M123(1,2)*B4 _b; 
Bl b=0; 
ifF>0 

X = -300e-10:deltax:0; 

Y = 0:deltax:L(1,2); 

Z = L(1,2):deltax:L(1,3); 

T = L(1,3):deltax:(L(1,3)+300e-10); 
else 

X = L(1,3):deltax:(L(1,3)+300e-10); 

Y = L(1,2):deltax:L(1,3); 

Z = 0:deltax:L(1,2); 

T = -300e-10:deltax:0; 
end 
P1_bound(a,:) = -C(1).*(F.*X + eta(1)); 
P2_bound(a,:) = -C(2).*(F.*Y + eta(2)); 
P3_bound(a,:) = -C(3).*(F.*Z + eta(3)); 
P4 _ bound(a,:) = -C(4).*(F.*T + eta(4)); 











psil_bound(a,:) = Al_b(:,a).*airy(0,P1_bound(a,:)) + 


B1 b.*airy(2,P1_bound(a,:)); 


B2_b(:,a).*airy(2,P2_bound(a,:)); 


psi2_bound(a,:) = A2_Db(:,a).*airy(0,P2_bound(a,:)) 1 





psi3_bound(a,:) = A3_b(:,a).*airy(0,P3_bound(a,:)) f 


B3_b(:,a).*airy(2,P3_bound(a,:)); 
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psi4_bound(a,:) = A4_b.*airy(0,P4_bound(a,:)) + B4_b.*airy(2,P4_bound(a,:)); 


psi_1_bound(a,:) = conj(psil_bound(a,:)); 
psi_2 bound(a,:) = conj(psi2_bound(a,:)); 
psi_3_bound(a,:) = conj(psi3_bound(a,:)); 
psi_ 4 bound(a,:) = conj(psi4_bound(a,:)); 


integrand! bound(a,:) = (psil_bound(a,:).*psi_1_bound(a,:))*deltax; 
integrand2 bound(a,:) = (psi2_bound(a,:).*psi_2 bound(a,:))*deltax; 
integrand3_bound(a,:) = (psi3_bound(a,:).*psi_3_bound(a,:))*deltax; 
integrand4 bound(a,:) = (psi4_bound(a,:).*psi_4 bound(a,:))*deltax; 


integrall_bound(a,:) = sum(integrand1 bound(a,:)); 
integral2_bound(a,:) = sum(integrand2_bound(a,:)); 
integral3_bound(a,:) = sum(integrand3_bound(a,:)); 
integral4 _bound(a,:) = sum(integrand4 bound(a,:)); 


finali_bound(a,:) = integrall bound(a,:) + integral2 bound(a,:) +  inte- 
gral3_bound(a,:) + integral4_bound(a,:); 


b3_bound(a,:) = 1/(sqrt(finali_ bound(a,:))); 


Psil_bound(a,:) = b3_bound(a,:)*psil_bound(a,:); 
Psi2_bound(a,:) = b3_ bound(a,:)*psi2_bound(a,:); 
Psi3_bound(a,:) = b3_ bound(a,:)*psi3_ bound(a,:); 
Psi4_bound(a,:) = b3_bound(a,:)*psi4_bound(a,:); 


%**Simulation Of Performance Of Infrared Photodetectors******* 
%Epsi_cont 

%LT Psarakis Eftychios Hellenic Navy 

deltax = 10-11; 


if F>0 
X = -300e-10:deltax:0; 
Y = 0:deltax:L(1,2); 
Z = L(1,2):deltax:L(1,3); 
T = L(1,3):deltax:(L(1,3)+300e-10); 


AA(ac) = 1/sqrt(2*pi); 

B4(ac) = -(M123_c(2,1)/M123_c(2,2))*A4(ac) ; 
A3(ac) = (M_c(1,3)*A4(ac) + B4(ac)*M_c(2,3)); 
B3(ac) = (M_c(3,3)*A4(ac) + B4(ac)*M_c(4,3)); 
A2(ac) = (M23_c(1,1)*A4(ac)+B4(ac)*M23_c(1,2)); 





o2 


B2(ac) = (M23_c(2,1)*A4(ac)+B4(ac)*M23_c(2,2)); 
Al(ac) = (M123 _c(1,1)*A4(ac) + B4(ac)*M123_c(1,2)); 
Bl(ac) = 0; 

else 


X = L(1,3):deltax:(L(1,3)+300e-10); 
Y = L(1,2):deltax:L(1,3); 

Z = 0:deltax:L(1,2); 

T = -300e-10:deltax:0; 


AA(ac) = 1/sqrt(2*pi); 
BA(ac) = -(M123_c(2,1)/M123_c(2,2))*A4(ac) ; 
A3(ac) = (M_c(1,3)*A4(ac) + B4(ac)*M_c(2,3)); 
B3(ac) = (M_c(3,3)*A4(ac) + B4(ac)*M_c(4,3)); 
A2(ac) = (M23_c(1,1)*A4(ac)+B4(ac)*M23_c(1,2)); 
B2(ac) = (M23_c(2,1)*A4(ac)+B4(ac)*M23_c(2,2)); 
Al(ac) = (M123 _c(1,1)*A4(ac) + B4(ac)*M123_c(1,2)); 
Bl(ac) = 0; 

end 








P1_cont(ac,:) = -C(1).*(F.*X + eta(1)); 
P2_cont(ac,:) = -C(2).*(F.*Y + eta(2)); 
P3_cont(ac,:) = -C(3).*(F.*Z + eta(3)); 
P4 cont(ac,:) = -C(4).*(F.*T + eta(4)); 








psil_cont(ac,:) = Al(ac)*airy(0,P1_cont(ac,:)) + Bl(ac)*airy(2,P1_cont(ac,:)); 
psi2_cont(ac,:) = A2(ac)*airy(0,P2_cont(ac,:)) + B2(ac)*airy(2,P2_cont(ac,:)); 
psi3_ cont(ac,:) = A3(ac)*airy(0,P3_cont(ac,:)) + B3(ac)*airy(2,P3_cont(ac,:)); 
psi4_cont(ac,:) = A4(ac)*airy(0,P4_cont(ac,:)) + B4(ac)*airy(2,P4_cont(ac,:)); 





psi_1_cont(ac,:) = conj(psil_cont(ac,:)); 
psi_2 cont(ac,:) = conj(psi2_cont(ac,:)); 
psi_3_cont(ac,:) = conj(psi3_cont(ac,:)); 
psi_4 cont(ac,:) = conj(psi4_cont(ac,:)); 


integrand! cont(ac,:) = (psil_cont(ac,:).*psi_1_cont(ac,:))*deltax; 
integrand2_ cont(ac,:) = (psi2_cont(ac,:).*psi_2 cont(ac,:))*deltax; 
integrand3_cont(ac,:) = (psi3_cont(ac,:).*psi_3_cont(ac,:))*deltax; 
integrand4 cont(ac,:) =(psi4_cont(ac,:).*psi_4 cont(ac,:))*deltax; 


integrall_cont(ac,:) = sum(integrand1_cont(ac,:)); 
integral2_cont(ac,:) = sum(integrand2_cont(ac,:)); 
integral3_cont(ac,:) = sum(integrand3_cont(ac,:)); 
integral4_cont(ac,:) = sum(integrand4 cont(ac,:)); 
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finali_cont(ac,:) =  integrall cont(ac,:) + integral2 cont(ac,:) +  inte- 
gral3_cont(ac,:) + integral4_cont(ac,:); 


b3_cont(ac,:) = 1/(sqrt(finali_cont(ac,:))); 

Psil_cont(ac,:) = b3_cont(ac,:)*psil_cont(ac,:); 
Psi2_cont(ac,:) = b3_cont(ac,:)*psi2_cont(ac,:); 
Psi3_cont(ac,:) = b3_cont(ac,:)*psi3_cont(ac,:); 
Psi4_cont(ac,:) = b3_cont(ac,:)*psi4_cont(ac,:); 


psil_cont_der(ac,:) = Al(ac)*airy(1,P1_cont(ac,:)) 
+Bl(ac)*airy(3,P1_cont(ac,:)); 

psi2_cont_der(ac,:) = A2(ac)*airy(1,P2_cont(ac,:)) + 
B2(ac)*airy(3,P2_cont(ac,:)); 

psi3_cont_der(ac,:) = A3(ac)*airy(1,P3_cont(ac,:)) + 
B3(ac)*airy(3,P3_cont(ac,:)); 

psi4_cont_der(ac,:) = A4(ac)*airy(1,P4_cont(ac,:)) + 
B4(ac)*airy(3,P4_cont(ac,:)); 

Psil_cont_der(ac,:) = b3_cont(ac,:)*psil_cont_der(ac,:); 

Psi2_cont_der(ac,:) = b3_cont(ac,:)*psi2_cont_der(ac,:); 

Psi3_cont_der(ac,:) = b3_cont(ac,:)*psi3_cont_der(ac,:); 

Psi4_cont_der(ac,:) = b3_cont(ac,:)*psi4_cont_der(ac,:); 

Integrall_cont(ac,:)}=sum((Psil_cont(ac,:).*conj(Psil_cont(ac,:)))*deltax); 

Integral2_cont(ac,:)=sum((Psi2_cont(ac,:).*conj(Psi2_cont(ac,:)))*deltax); 

Integral3_cont(ac,:)=sum((Psi3_cont(ac,:).*conj(Psi3_cont(ac,:)))*deltax); 

Integral4_cont(ac,:)}=sum((Psi4_cont(ac,:).*conj(Psi4_cont(ac,:)))*deltax); 


To- 
tal_integral_cont(ac,:)=Integrall_cont(ac,:)+Integral2_cont(ac,:)+Integral3_cont(ac,:)+Int 
egral4_cont(ac,:); 


%**Simulation Of Performance Of Infrared Photodetectors******* 
%LT Psarakis Eftychios Hellenic Navy 

function y=aproxAiry(x) 

“%othis function aproximates the airy's functions and their 
“odirivatives 

% Author: © E.Psarakis, § 31-Aug-2004 
%Ai(x)-->ans(1,1) --> ai 

%Bi(x)-->ans(1,2) --> bi 

%Ai'(x)-->ans(1,3)--> ad 

%Bi'(x)-->ans(1,4)--> bd 

Yoai=[];bi=[];ad=[];bd=[]; 
%ai=(4*pi)(-1/2)*(x).4(-1/4).*exp(-2/3*(x).%(3/2)); 
Yobi=(p1)(-1/2)*(x).(-1/4).*exp(2/3*(x).4(3/2)); 
%ad=-(4*pi)*(-1/2)*(x).(1/4). *exp(-2/3 *(x).4(3/2)); 
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%bd=(pi)(-1/2)*(x).4(1/4).*exp(2/3*(x).4(3/2)); 
%[ai; bi sad ;bd]; 
“%y=eval('[ai;bi;ad;bd]'); 
index=0; 
ai=[];bi=[];ad=[];bd=[]; 
for index=1 :length(x); 

if x(index)>=0 
ai(index)=(4*pi)(-1/2)*(x(index)).(-1/4).*exp(-2/3*(x(index)).“1.5); 
bi(index)=(p1)(-1/2)*(x(index)).(-1/4).*exp(2/3*(abs(x(index))).%1.5); 
ad(index)=-(4*pi)*(-1/2)*(x(index)).%(1/4).*exp(-2/3*(x(index)).%1.5); 
bd(index)=(pi)*(-1/2)*(x(index)).4(1/4).*exp(2/3*(x(index)).“1.5); 
else 
ai(index)=(pi)*(-1/2)*(-(x(index))).“(-1/4).*sin(2/3*(-(x(index))).%1.5+p1/4); 
bi(index)=(p1)(-1/2)*(-(x(index))).(-1/4).*cos(2/3*(-(x(index))).1.5+pi/4); 
ad(index)=-(pi)(-1/2)*(-(x(index))).4(1/4).*(cos(2/3 *(-(x(index))).1.5+p1/4)); 
%-1/(4*sqrt(pi).*(abs(x(index))).%(-5/4)).*sin(2/3*(abs(x(index))).%1.5+p1/4); 








bd(index)=(pi)(-1/2)*(-(x(index))).4(1/4).*(sin(2/3*(-(x(index))).“1.5+pi/4)); 
%-1/(4*sqrt(pi).*(abs(x(index))).“(-5/4)).*cos(2/3*(abs(x(index))).*1.5+pi/4); 
end 

end 

[ai; bi;ad ;bd]; 

y=eval('[ai; bi;ad ;bd]'); 

% ai=[];bi=[];ad=[];bd=[]; 








%**Simulation Of Performance Of Infrared Photodetectors******* 
%LT Psarakis Eftychios Hellenic Navy 
function t=asymptoticAiry(n,x) 
Y%othis function aproximates the airy's functions and their 
“%dirivatives. it is used like the matlab airy(x) 
% Author: © E.Psarakis, § 31-Aug-2004 
%Ai(x)-->ans(1,1) 
%Bi(x)-->ans(1,2) 
%Ai'(x)-->ans(1,3) 
%Bi'(x)-->ans(1,4) 
if n==0 
aprox Airy(x); 
t=ans(1,:); 
end 
if n== 
aprox Airy(x); 
t=ans(3,:); 
end 
if n== 
aprox Airy(x); 
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t=ans(2,:); 
end 
if n==3 

aprox Airy(x); 

t=ans(4,:); 
end 


%**Simulation Of Performance Of Infrared Photodetectors******* 
%ZETA 

%LT Psarakis Eftychios Hellenic Navy 

warning off 

deltaz = 10-11; 


ifF >0 

X = -300e-10:deltaz:0; 

Y = 0:deltaz:L(1,2); 

Z = L(1,2):deltaz:L(1,3); 

T = L(1,3):deltaz:(L(1,3)+300e-10); 

dim=max(T)-min(X); 
else 

X = L(1,3):deltaz:(L(1,3)+300e-10); 

Y = L(,2):deltaz:L(1,3); 

Z = 0:deltaz:L(1,2); 

T = -300e-10:deltaz:0; 

dim=max(X)-min(T); 
end 
for kk=1:ac-1 
Zital(kk,:)=conj(Psil_bound(1,:)).*X.*Psil_cont(kk,:).*deltaz; 
Zita2(kk,:)=conj(Psi2_bound(1,:)).*Y.*Psi2_cont(kk,:).*deltaz; 
Zita3(kk,:)=conj(Psi3_bound(1,:)).*Z.*Psi3_cont(kk,:).*deltaz; 
Zita4(kk,:)=conj(Psi4_bound(1,:)).*T.*Psi4_cont(kk,:).*deltaz; 
OSC_integral1(kk,:) = sum(Zital(kk,:)); 
OSC_integral2(kk,:) = sum(Zita2(kk,:)); 
OSC_integral3(kk,:) = sum(Zita3(kk,:)); 
OSC_integral4(kk,:) = sum(Zita4(kk,:)); 
OSC_finali(kk,:)}=(OSC_ integral 1(kk,:)+ OSC_integral2(kk,:)+ 

OSC_integral3(kk,:)+OSC_integral4(kk,:)); 

Zeta_final(kk,:)=(OSC_finali(kk,:)).*conj(OSC_finali(kk,:)); 


% DOS(kk)=dim./(hbar*pi)*sqrt(m_eff(1)./(2*q*(Eout_cont(kk)- 
min(V(1),V(2))))); 

DOS(kk)=dim./(hbar*pi)*sqrt(m_eff(1)./(2*q*(Eout_cont(kk)- 
Eout_cont(1)+.001))); 


end 
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APPENDIX B: EKF INITIALIZATION AND TRACKING 


In the present Appendix we present the BOT initialization and tracking algorithm. 
The core of the algorithm is BOT.m. Once this program is executed, the loop INITIALI- 


ZATION.m is called, and the user must define the initialization bearings. 


MW***EEEWriten by LT Eftychios Psarakis Hellenic Navy*******#% #2 2s ok ok 


%***Thesis Project: Simulation of Performance of Infrared Photodetectors***** 
1/2 Fe 6 He 2 2 2 ae he 8 2 2 2 fe ie 2 2 2 6 6 he 2 2 ae fe ie 2 2 2 ae 46 Ae 2 2 ee 6 26 2 2g ae fe 6 2 2 2 fe ie 26 2 2s ae fe ie 2 2s 2c fe ie 2s 2s 2s ae 


% cle 

% clear 

Slp=[0;0]; % Sensor | at origin 
S1p=[10000;0]; % Sensor 2 at (10000,0) 


% bl1=input('give mesured bearing | from sensor | in degrees: ');%33 
% b2=input('give mesured bearing | from sensor 2 in degrees: ');%173 
bl_1=b1*2*pi/360; Yoconvert to rad 

bl_2=b2*2*pi/360; Yoconvert to rad 

sigma 1b=1*pi/180; % 1 deg std dev in bearing 
miss=abs(randn*sigma1b); 


b1_S1_plus =b1_1+miss; % bearing 1 from sensor 1 (+) 
bl_S1_minus=b1_1-miss; % bearing | from sensor 1 (-) 


bl_S2_plus =b1_2+miss; % bearing | from sensor 2 (+) 
b1_S2_minus=b1_2-miss; % bearing 1 from sensor 2 (-) 


% b_1=input('give mesured bearing 2 from sensor | in degrees: ');%33.46 
% b_2=input('give mesured bearing 2 from sensor 2 in degrees: ');%173.18 
% dt=input('give time interval between the two measurements : '); 

dt=.1; 

b2_1=b_1*2*pi/360; Yoconvert to rad 

b2_2=b_2*2*pi/360; Yoconvert to rad 


b2_S1_plus=b2_1+miss; % bearing 2 from sensor | (+) 
b2_S1_minus=b2_1-miss; % bearing 2 from sensor | (-) 








b2_S2_ plus=b2_2+miss; % bearing 2 from sensor 2 (+) 


b2_S2_minus=b2_2-miss; % bearing 2 from sensor 2 (-) 
% Ce ee 


%Triangle 1: Bearing | Sensor 1(+), Bearing 1 Sensor 2(+) using sine law 
theta3_1=pi-b1_S1_plus-(pi-b1_S2_ plus); 


a7 


D1=((sin(theta3_1)/10000)*-1)*sin(pi-b1_S2_ plus); Yodistance from (0,0) 


xl=D1*cos(b1_S1_ plus); 
yl=D1*sin(b1_S1_ plus); 


% Triangle 1: Bearing 1 Sensor 1(+), Bearing 1 Sensor 2(-) using sine law 
theta3_2=pi-b1_S1_plus-(pi-b1_S2_ minus); 
D2=((sin(theta3_2)/10000)*-1)*sin(pi-b1_S2_ minus); %distance from (0,0) 
x2=D2*cos(b1_S1_ plus); 

y2=D2*sin(b1_S1_plus); 


% Triangle 1: Bearing 1 Sensor 1(-), Bearing 1 Sensor 2(+) using sine law 
theta3_3=pi-b1_S1_minus-(pi-b1_S2_ plus); 
D3=((sin(theta3_3)/10000)*-1)*sin(pi-b1_S2_ plus); Yodistance from (0,0) 
x3=D3*cos(bl_S1_minus); 

y3=D3*sin(b1_S1_ minus); 


%Triangle 1: Bearing 1 Sensor 1(-), Bearing 1 Sensor 2(-) using sine law 
theta3_4=pi-b1_S1_minus-(pi-b1_S2_minus); 
D4=((sin(theta3_4)/10000)*-1)*sin(pi-b1_S2_ minus); %distance from (0,0) 
x4=D4*cos(b1_S1_ minus); 

y4=D4*sin(b1_S1_ minus); 


% disp('*****#*** first parallelogram ********##*+44+!): 
posl=[xl;y1]; 
pos2=[x2;y2]; 
pos3=[x3;y3]; 
pos4=[x4;y4]; 
% displ tee eee ne en ee Tee ae) 
POS 1=[(x1+x2+x3+x4)/4; 
(yl+y2+y3+y4)/4]; % middle of the rectungular shape 


> 


d=max(abs((x2-x3)/(atan((y2-y3)/(x2-x3)))),abs((x4-x1)/(atan((y4-y1)/(x4- 
x1))))); 


xx=linspace(1,10000,100); 
yyl=tan(b1_S1_plus)*xx; 
yy2=tan(b1_S1_minus)*xx; 


a3=y3/(x3-10000); 
b3=-10000*a3; 
yy3=a3*xx+b3; 
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a4=y4/(x4-10000); 
b4=-10000*a4; 
yy4=a4*xx+b4; 
% figure 

% plot(xx,yy1,'r') 
% hold on 

% plot(xx,yy2,'b') 
% hold on 

% plot(xx,yy3,'r') 
% hold on 

% plot(xx,yy4,'b') 


% Triangle 2: Bearing 1 Sensor 1(+), Bearing 1 Sensor 2(+) using sine law 
theta3_1_2=pi-b2_S1_plus-(pi-b2_S2_ plus); 
D1_2=((sin(theta3_1_2)/10000)*-1)*sin(pi-b2_S2_ plus); %distance from (0,0) 
xl_2=D1_2*cos(b2_S1_ plus); 

yl_2=D1_2*sin(b2_S1_ plus); 


% Triangle 2: Bearing 1 Sensor 1(+), Bearing 1 Sensor 2(-) using sine law 
theta3_ 2 2=pi-b2_S1_plus-(pi-b2_S2_minus); 

D2_2=((sin(theta3_2 2)/10000)*-1)*sin(pi-b2_S2 minus); %odistance from (0,0) 
x2_2=D2_2*cos(b2_S1_plus); 

y2_2=D2_2*sin(b2_S1_ plus); 


% Triangle 2: Bearing 1 Sensor 1(-), Bearing 1 Sensor 2(+) using sine law 
theta3_3_2=pi-b2_S1_minus-(pi-b2_S2_plus); 
D3_2=((sin(theta3_3_2)/10000)*-1)*sin(pi-b2_S2_ plus); “distance from (0,0) 
x3_2=D3_2*cos(b2_S1_ minus); 

y3_2=D3_2*sin(b2_S1_ minus); 


% Triangle 2: Bearing 1 Sensor 1(-), Bearing 1 Sensor 2(-) using sine law 
theta3_4 2=pi-b2_S1_minus-(pi-b2_S2_ minus); 

D4_2=((sin(theta3_ 4 2)/10000)*-1)*sin(pi-b2_S2 minus); %odistance from (0,0) 
x4 2=D4 2*cos(b2_S1 minus); 

y4_2=D4_2*sin(b2_S1_ minus); 


% disp('*********second parallelogram” ******* #%## x #71). 
posl 2=[x1_2;yl_2]; 
pos2_2=[x2_2;y2_2]; 
pos3_2=[x3_2;y3_ 2]; 
pos4 2=[x4 2;y4 2]; 


% igp(' tPF ttt Eee eee LST ELE ESE SESE ES SEERA EEE EEE: 


POS2=[(x1_2+x2_2+x3_2+x4 2)/4; 
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(yl _2+y2_ 2+y3 2+y4 2)/4]; % middle of the rectungular shape 
dl=max(abs((x2_2-x3_2)/(atan((y2_2-y3_2)/(x2_2-x3_2)))),abs((x4_2- 
xl _2)/(atan((y4_2-yl_2)/(x4_2-x1_ 2))))); 











sigma 1r=max(d,d1); 
xx=linspace(1,10000,100); 


yyl_2=tan(b2_S1_plus)*xx; 
yy2_2=tan(b2_S1_minus)*xx; 


a3_2=y3_2/(x3_2-10000); 
b3_2=-10000*a3_2; 
yy3_2=a3_ 2*xx+b3_2; 


a4_2=y4 2/(x4_2-10000); 
b4_2=-10000*a4 2; 
yy4_ 2=a4 2*xxt+b4 2; 
% figure 

% plot(xx,yy1_2,'r') 

% hold on 

% plot(xx,yy2_2,'b') 

% hold on 

% plot(xx,yy3_2,'r') 

% hold on 

% plot(xx,yy4_2,'b') 


x1=[POS2(1,1);POS2(2,1);POS1(1,1);POS1(2,1)]; 
% 

% figure 

% plot(xx,yy1,'r') 

% hold on 

% plot(xx,yy2,'r') 

% hold on 

% plot(xx,yy3,'r') 

% hold on 

% plot(xx,yy4,'r') 

% hold on 

% plot(xx,yy1_2,'b') 
% hold on 

% plot(xx,yy2_2,'b') 
% hold on 

% plot(xx,yy3_2,'b') 
% hold on 
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% plot(xx,yy4_2,'b') 

% 

% text(POS1(1,1),POS1(2,1),'POS1') 
% text(POS2(1,1),POS2(2,1),'POS2') 


Oo) 7 ae He 2k 2 2s ae he he 2 2g 2g fe 6 26 2 2 2S AEB YT yy 2 2 2 2g ie ie 2 2 2g he fe he 2 2 2 6 he 26 2 2g 2g fe ie 2 2 2g ae fe ie 2 2s 2s Ie ik 


Wo sR riten by LT Eftychios Psarakis Hellenic Navy********** 
%***Thesis Project: Simulation of Performance of Infrared Photodetectors*** 
1/7 2s 2s ae he he 2 He 2g fe hee 2 2 ae 6 6 2 2 2 6 6 26 2 2 6 fe 6 2 2 2 6 46 Ae 2 2 2 fe 6 26 2 2 6 6 Ae 2 2s 2g fe ie 2 2 2s ie fe ie 2 2s 2s fe ie 
clear,clc 

disp(‘for no mesurement noise select 0') 

disp('for mesurement noise select 1") 

mnoiseflag=input(“SELECT NOISE OR NO NOISE CASE::); 


disp(‘for single sensor tracking select 0') 
disp(‘for two sensor tracking select 1') 
twosensorflag=input('SELECT ONE OR TWO SENSORS CASE::’); 


s2p =[10000;0];  % Sensor 2 position (Sensor | is at (0/0) 
%rand('seed',0); , 

delta = 0.1; %Time 

nsamples = 1000; 


initialization; 
X =x 
q = 1000; 


Q= [delta*3/3, delta*2/2; delta*2/2, delta]; 
Qi = [Q, zeros(2); zeros(2),Q]; 


Q= q*Qi; 

sb = 1*pi/180; 

D=[1, 0, 0, 0; 
1/dt, 0, -1/dt, 0; 
0, 1, 0, 0; 
0, 1/dt, 0, -1/dt]; 

x_sp=D*x1; 


speed=sqrt((x_sp(2,1))42+(x_sp(4,1))*2)*2000/3600 
w1=input('SELECT THE DESIRED ACCELERATION IN G') 
w = wl*10/speed; 


R=sb; 

xout = []; 

tout = []; 

terr = []; 

chi2err = []; 

% xhat = x + [100*randn;5*randn;100*randn;5*randn]; 
xhat = x + 0.6*[100;5;100;5]; 


101 


phat=[d1’2 ,0 ,0 ,0; 


0 ,d1%2,0_ ,0; 
O.~ 30. gd*2.50: 
0 0 0 ,d%2] 
F=[0, I, 0, 0; 
0, 0, 0 -W; 
0, 0, 0, 1; 
0, WwW, 0, 0]; 
Fturn = expm(F*de1ta) 
F =[1, delta, 0, 0; 
0, i 0 0; 
0, 0, 1, delta; 
0, 0, 0, 1); 


Hp =[1000;00 1 0]; 
for 11 =1:nsamp les 
% Take a measurement 


b = atan2(x(3),x(1)); 

if mnoiseflag== 
z=b+sb*randn; 

else 

z=b; 

end 

xout = [xout, [x(1) ;x(3) ] ]; 


% EKF Update 

r2 =xhat(1)42 + xhat(3)%2; 

b = atan2(xhat(3),xhat(1)); 

zhat = b; 

H= [-xhat(3)/r2, 0, xhat(1)/r2,0]; 
Pzi=inv(H*phat*H' + R); 

zt = z - zhat; 

chi2 =zt'*Pzi*zt; 

K = phat*H'*Pzi; 

K2 = eye(4)-K*H; 

xhat =xhat + K*zt; 

phat = K2*phat*K2' + K*R*K’; 


if twosensorflag ==1 % process measuransnt frcm second sensor 
% Take a Measuronent 
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rp =Hp*x - s2p; 

b = atan2(rp(2),rp(1)); 

if mnoiseflag == 
z=b+sb*randn; 

else 

zZz=b; 

end 

“% Measuranent update 

rp = Hp*xhat - s2p; 

12 = rp'*rp; 

b = atan2(rp(2),rp(1)); 

zhat = b; 

H= [-rp(2)/r2, 0,rp(1)/r2,0] ; 
Pzi = inv(H*phat*H' + R) ; 
zt = z- zhat; 

chi2 = zt'*Pzi*zt; 

K = phat*H'*Pzi; 

K2 = eye(4) - K*H; 

xhat = xhat + K*zt; 

phat = K2*phat*K2' + K*R*K’; 


end 

e = [xhat(1);xhat(3)]; 

tout =[tout,e] ; 

e=e - [x(1);x(3)]; 

terr = [terr,sqrt(e'*e)]; 

chi2err = [chi2err,chi2]; 

% Target Motion 

if (1 > 250) & (ii <= 350)) 

x = Fturn*x; 

else 

x = F*x; 

end 

% Track Prediction 

xhat=F *xhat; 

phat=F*phat*F' + Q; 

end % for 11 

if mnoiseflag==0 && twosensorflag==0 

tlabel = ['BO Tracking single Sensor,no Noise,q‘2 = 1000, accel 
='~num?2str(w1),'g"]; 

elseif mnoiseflag==1 && twosensorflag==0 

tlabel = ['BO Tracking single Sensor,Noise presence,q*’2 = 1000, accel 
='"num2?2str(w1),'g']; 

elseif mnoiseflag==0 && twosensorflag== 

tlabel = ['BO Tracking two Sensors, no Noise,q®2 = 1000, accel 
='~num?2str(w1),'g"]; 
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elseif mnoiseflag==1 && twosensorflag== 

tlabel = ['BO Tracking two Sensors,Noise presence,q*2 = 1000, accel 
='"num2?2str(w1),'g']; 

end 

figure 

subplot(2,1,1) 

plot(xout(1,:),xout(2,:), ':',tout(1, :),tout(2,:),'r'); 

%title( 'Target and Track’) ; 

title ([tlabe1]) ; 

xlabel('meters') 

ylabel('meters') 

legend('true position','Track position’) 

grid 

subplot(2,1,2) 

plot(xout(1,:),xout(2,:), ",tout(1, :),tout(2,:),'r'); 

%title( 'Target and Track’) ; 

title ([tlabe1]) ; 

xlabel('meters') 

ylabel('meters') 

legend('true position','Track position’) 

axis([.6*10%4,3.3*10%4,3*10%4,13*10%4]) 

grid 

figure 

tt = [0: (max (size(terr))-1)] *del ta; 

subplot(2,1,1) 

plot (tt,terr, 'r'); 

title ({ 'Track Position Errors: ',tlabe1]) ; 

xlabel('Sample number") 

ylabel('meters') 

grid 

subplot(2,1,2) 

tt = [0: (max(size (chi2err))-1)] *delta; 

plot (tt,chi2err, 'r') ; 

title (‘Chi-Squared Values for Measurement Association’); 

xlabel('Sample number") 

ylabel('rad"') 

grid 


104 


APPENDIX C: ASYMPTOTIC EXPRESSIONS 


Appendix C presents the Airy function considerations for the low bias simulation 
of the quantum well. 
A. AIRY FUNCTION ANALYSIS 

1. General 


2 
The differential equation > v0 = 0, where  =2/3z° andw=e””, can be 
Iz 


reduced to the differential equation satisfied by Bessel functions of order 1/3 [33, sec. 
6.4]. There are two linearly independent solutions which are called Airy functions of the 
first and the second kind. 


(j= Ke [as($)- his (6) (C.1) 
! AGL to Looleehelé (C.2) 
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Figure C.1: Airy functions of the first kind, Ai(—z) and Ai(z) 
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Airy functions of the second kind 














Figure C.2: Airy functions of the second kind, Bi(—z) and Bi(z) 


oe Asymptotic Expressions of Airy Functions 


There are several problems associated with the direct approach to the solution of 
the Schrédinger equation using Airy’s functions. The main drawback is the asymptotic 
behavior of Airy’s functions and of their derivatives, when the used argument, is either 
very large or very small. In those cases personal computers give numerical overflows 
which significantly reduce the result quality and computational efficiency. These kinds of 
problems appear especially when low biases are applied in symmetric or asymmetric 
quantum wells structures. 

The asymptotic forms of Airy’s functions used for our simulations were taken 


from the handbook of mathematical functions [28, pp. 448-449]. 


a ae: 
Ai(z)= (4m) 2z 4e? ,z>>0 (C.5) 
Da 25 
Ai (z) = -(42) 2z4e 3 ,z>>0 (C.6) 
& ft  f5. 2 g 
Ai(z) = (2) 7(—z) 4 snl 22) 2) << 0 (C.7) 
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Ai (z) = -(2) ?(-z)4 cos 2-0) + z z<<0 (C.8) 


HORE tLe eS (C.9) 
Bi(z)=(a)2 er »2380 (C.10) 
Bi(z) = (2)? eae cs( 2-0) + z ,z<<0 (C.11) 
Bit) =(a) *Catsin[ 2-0) +2) 20 (C.12) 


The expressions above give acceptable results not only when the argument is very 
large or extremely small but also when the argument is a finite number (not around zero) 
as shown in Figures C.3 and C.4. In the simulation code, written in Matlab, we have used 
a specific procedure so that when the built-in Matlab Airy functions give results, we take 
advantage and use them, but when the built-in functions are unable to give results we use 


Equations (C.5) through (C.12). 


Matlab Airy Vs Asymptotic Expressions of the first kind 
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Figure C.3: Matlab Airy functions vs. asymptotic expressions of the first kind 
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Figure C.4: Matlab Airy functions vs. asymptotic expressions of the second kind 
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BOT COMPLETE CASE SET 


APPENDIX D 


We present the complete results of the BOT simulation for 0 g, 2 g, 4 g, 6 g and 8 


g turns, using one or two sensors, with and without the presence of measurement noise. 
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